{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pylab as pl\n",
    "import itertools\n",
    "\n",
    "pl.rcParams[\"figure.figsize\"] = 4, 4\n",
    "\n",
    "from gmm_base import *\n",
    "from simulation import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "background = 0.2\n",
    "base = 2\n",
    "lamb = lambda t: background + base**-t\n",
    "cond = 1\n",
    "xmax = 3\n",
    "bins = np.linspace(0, 10, 201)\n",
    "sample_size = int(1e4)\n",
    "num_clusters = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### data generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array(\n",
    "    [\n",
    "        thinning_sampler(rng, lamb)\n",
    "        for rng in itertools.repeat(np.random.RandomState(0), sample_size)\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 3)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAD4CAYAAAAZ+NgoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2dd3hUVfrHP2cmE9IhJLQUEpokkEZHelOQdQEREJS2otjXVVHx565iR2VddcUOAoqiYkOWIlXpHUxC0dBDFAi9pE7O74+bDJNkktwkM5mS83me+zjnnnPPfe819+WU93yPkFKiUCgUejA42wCFQuE+KIehUCh0oxyGQqHQjXIYCoVCN8phKBQK3Xg568ahoaEyOjraWbdXKGotO3bsyJRSNqjKtU5zGNHR0Wzfvt1Zt1coai1CiKNVvVZ1SRQKhW6Uw1AoFLpRDkOhUOjGaWMYnkZeXh7p6elkZ2c72xSFAgAfHx8iIiIwmUx2q1M5DDuRnp5OYGAg0dHRCCGcbY6iliOl5MyZM6Snp9OsWTO71au6JHYiOzubkJAQ5SwULoEQgpCQELu3eCt0GEKI2UKIU0KIlDLy7xBC/Fp4bBRCJNrVQjdCOQuFK+GIv0c9LYw5wKBy8g8DvaWUCcALwId2sEuhULggFToMKeUvwNly8jdKKc8VJjcDEXpufCErD6XFoVC4F/Yew5gELC0rUwgxWQixXQix/djZqxw/m2Xn2ysUjufcuXMVF/JQ7OYwhBB90RzGk2WVkVJ+KKXsKKXs2KphAGH1fOx1e4WixnjkkUecbYLTsIvDEEIkAB8DQ6WUZ/Rc42My4mVUkzQ1RVZWFr1798ZsNgPaNPCXX35Z7jXJyclERUXx3nvv6a4rNzeXXr16kZ+fX2Vbly1bRuvWrWnZsiXTp08vlX/8+HH69u1LbGwsbdu25a233qryvSpb17Jly9i/fz8zZsywy7O6G9X+YoUQTYFvgXFSyt/0Xnc118wXW49V9/YKncyePZvhw4djNBoBWLVqFTt37iz3mvj4eBYsWMC8efN01+Xt7U3//v0rdEZlYTabeeCBB1i6dCl79+7liy++YO/evcXKeHl58e9//5t9+/axefNmZs6cWaqMNWvXrmXixIk288qrKzk5mZtvvrnYUb9+fcaOHcuUKVOq/axuiZSy3AP4AvgDyAPS0bod9wL3FuZ/DJwDdhce2yuqU0pJWMu2MuafS2Vevll6Anv37nW2CVJKKW+77TY5atQo2blzZ9m0aVO5ePFiKaWU119/vTx8+LCUUsp169bJ+vXry+bNm8vExER56NAhedNNN8kTJ06Uqi8tLU0GBgYWO1dRXbt375Y33XRTlezfuHGjvPHGGy3pl19+Wb788svlXjNkyBD5008/lZm/Zs0aOWHCBF33r6iuWbNmyV9++cWSrs6z1gS2/i71fqO2jipdZI8jsV17eSk7z06vxfm4isOIiYmRU6dOlVJqH3OnTp1kTk6ObNSoUbFyAwcOlMnJyRXWN2LECOnt7S2PHDkipZS66srPz5ehoaGl6urRo4dMTEwsdaxYscJS5uuvv5aTJk2ypOfNmycfeOCBMu07fPiwjIyMlBcuXCizjF6HoaeuH374QY4fP97y/7usZ3UV7O0wnBYafvhCGtcvaGdJh/mHsXzEcmeZY3du+2ATIzpEMLJjJHnmAsZ+vIXRnSO5pV0EWblmJn6ylbFdo/hrYhgXs/O4e+52/tY9mkFxTTh7JZf7PtvB3T2bM6BNI05dyqZhYMUDxFlZWWRmZvLss88C0KZNG86dO0dmZib16tUrVvbAgQO0bt263PqWLVvGlStX+Mtf/kJqaipRUVG66jIajXh7e3Pp0iUCAwMt59etW1fhM0gbU+1lBSBdvnyZW2+9lTfffJOgoKBS+V26dCEnJ4fLly9z9uxZkpKSAHj11VcZOHBgpeoqYsiQIQwZMsSSLutZPRWnOYzcglz+1WYJOeYCxnWNIn5uvLNM8RhSUlJo1aoVPj6ac9m5cyeJiYn4+voWCxE+c+YMdevWLXdRUnZ2Nk888QSLFi3ik08+ISUlhcGDB+uuKycnx2JHET179uTSpUul7jVjxgwGDBgAQEREBMePH7fkpaenExYWVuqavLw8br31Vu644w6GDx9u8xm2bNkCaGMYc+bMYc6cOTbL6amrPGw9q8dS1aZJdQ+faB85ac5WedsHG6WUUsbNibNDA8x5uEKX5KOPPpLh4eEyKytLXr58WXbr1k2uX79eSillRESEzMrKklJKuW3btlL97n79+sn09HRL+umnn5avv/66lFLrJowbN86SV1FdmZmZMiYmpkrPkJeXJ5s1ayYPHTokc3JyZEJCgkxJSSlWpqCgQI4bN04+/PDDuuosr0tS2bpKUp1nrQns3SVx6rzmO7e3Z8Hk651pgkexZ88e7rjjDvr06UOnTp2477776N69OwA33ngj69evByAmJobMzEzi4uLYuHEjBQUFpKWlUb9+fUDrYqxYsYJ//OMfgDZbkpJybSlReXUBrFmzhsGDB1fpGby8vHjnnXcYOHAgsbGxjBo1irZt2wIwePBgMjIy2LBhA59++imrV68mKSmJpKQklixZUqX7Vbeu6jyrW1JVT1Pdwyfap5jXUy2M6tOzZ0+5f/9+m3k7d+6UY8eOtZmXnJwsH3nkEd33Ka8uKaW85ZZbyrTD03D1Z/WYQU+AggLJtB9TiQ+v60wzPIaDBw/SqlUrm3nt2rWjb9++mM1mS/xEEXFxcbzxxhu671NeXbm5uQwbNqzCAVVPoDY9axFOdRgGg+DX9AvU9bWfIlBt5sSJE+Xm33nnnXa7V1l1eXt7M378eLvdx5WpTc9ahNMVt75/QOtjz5nrZEMUCkWFON1hFBHmH1ZsatXT4jIUCk/A6Q7j1MVsHvlqN1O6zeOGNo0s51VchkLhejh9uWg9P2+u5prJNxc42xSFQlEBTm9heHsZ+O7+7s42Q6FQ6MDpLYwiiuZ5FQqF6+ISDiPlxAW6TV/N5kNlSocqFAoXwCUcRmSwHx2igvHzNlZcWKFQOA2nj2EA1PUz8c7t7Z1thkLhUpw7d47g4GBnm1EMl2hhFHE5Jx9zgRrHUCjANcWGXcZhrDlwioRpy0nNuOBsUzySikSA+/Tpw5EjR2zmVZY777yThg0bEhcXp6u8Ehu2bZsrig27jMNo2ySIB/u1ItjP29mmeCSVEQEuK0/vfhwTJ05k2bJlum1TYsNuJDZc1WWu1T1KLm8vibstd3eF5e1SVl0EuHfv3vLw4cM284po3ry5HDNmjFy1apUsKCgo147Dhw/Ltm3bljqvxIb11WUvsWGPWt5eEnOBZP+fF2nTpGxNRbfhk7+UPtd2GHS+G3KvwvyRpfOTbod2d8CVM/BViVWQf/ufrtvu2bOHYcOG8eWXX7J+/XoeffRRbrjhBg4dOkR0dDQAPXr0oFOnTsyYMaNUt6G8vN9++42lS5fyzjvv8MADDzBu3DgmTpxoU0KvLMoSp5k6dSo5OTkcPXqUqKgocnNzK7TZbDazbdu2UnXpkQI8ceIEkZGRlryIiAiLpJ8tjhw5wq5du+jSpYvuZ61OXaGhoXz88ceEhoYSGxtLXFyczWetaVzKYSzccZwnv0lm1WO9nW2KW2IvEeCy8oxGo6XJfPr0aZ566imaNm3Kxo0b6dy5c5XtVmLDpXFVsWGXchi9r2vIW6OTaBBYx9mmVJ/yWgTefuXn+4foblFYYw8R4IoEgi9cuMCXX37JJ598gslkYtasWSQkJFTa1iKU2LB+XEFs2KUcRuO6PgxNCne2GW7Lnj17OHbsGNnZ2ZjNZp599llee+01goODMZvNZGdn4+Pjw+HDh8vsRpSXN3bsWDZt2sTIkSOZN29emepe5dG/f3/mzZtHeLj2//nFF19k/PjxREdHEx8fz6JFiwB02XzmzBkaNGhQyonoaWF06tSJ33//ncOHDxMeHs6CBQv4/PPPi5WRUjJp0iRiY2N59NFHK/2s9qyrrGetaVxmlqSI05dyWJ76p7PNcEuqKgJsTXl5o0aN4sCBA0yfPr1cZzFmzBiuv/56Dhw4QEREBLNmzQJQYsOeIDZc0agoMBs4BaSUkS+At4E04FegvZ7R1rJmST765aCMenKxbPNRt0qPCDsTV5glqaoIsJTSMkviSJTYcNWp6rM6Y5uBOcCgcvJvAloVHpOB98opWyFDEsP48cEeCK/L1ammVqJXBNhZVEdsuCS1SYDXlZ61wjEMKeUvQojocooMBeYVeq7NQoh6QogmUso/qmJQwyAfGgb5IIQS1Kks1REBnjhxYqlZCVdAiQ271rPaY9AzHDhulU4vPFfKYQghJqO1QjCFlD14s//Pi/hcutki06f0PR1PWRGKCoU19hj0tDV5bXMFmZTyQyllRyllR2Ng2UvZtxw6y5kTPVgxbBvJE5LJuJJhBzMVCkV1sUcLIx2ItEpHANX6wm9pH86wduFqvxKFwsWwRwtjETBeaHQFLlR1/KKIIB+TchYKhQtSocMQQnwBbAJaCyHShRCThBD3CiHuLSyyBDiENq36EXC/PQzbfOgMj3+9hwKlj6FQuAx6ZknGVJAvgQfsZlEhGeez+OX30/x5MbviwgqFokZwqdBwa4YkhnFLu/AyFwQpFIqax2UdhpfR5aLWFYpaj8s6DICfUv/kvZ8PIn2VmrhC4Qq4tMMwGQ34eBmRZn9nm1JpBi4caNf4EWcHr3Xr1q3UYjQ9eeUREBDA5cuVWwJw/vx5Pv/8c+6/X9/YepH83dSpU7nvvvss57Oyshg0aBCrV69mxYoVPPzww5jNZu666y6mTp0KaCHZAwYMYPXq1Xh5Ve1TufPOO1m8eDENGzYstrDOmuPHjzN+/Hj+/PNPDAYDkydP5uGHHwY0rRBbtpV13uFUdRFKdY+KJPqscQe5vpKLfOxtsyu+g4KCAmk2m6t8vb+/f6WvKUv6rzw2btwou3btWuzcO++8I998802Zn58vmzdvLg8ePChzcnJkQkKCTE1NtZSbNm2a/OyzzyptZxE///yz3LFjR7k2Z2RkyB07dkgppbx48aJs1aqVTE1NLdO2imy2xhmLz5yOlGrgUy+fffYZnTt3JikpiXvuuQez2cyRI0eIiYnhrrvuIi4ujjvuuIOVK1fSvXt3WrVqxdatWy1lJkyYQEJCAiNGjODq1auWegMCAgBNXi42Npb777+f9u3bc/z4cUvevHnzSEhIIDExkXHjxlmuHTZsGB06dKBt27Z8+OGH5dpfpLyVmJhIXFxcKeHbqVOncvDgQZKSknj88cd1vZOGDRuSmppa7Nz8+fMZOnQoW7dupWXLljRv3hxvb29Gjx7NDz/8UMz2+fPn67qPLXr16mVZzl8WTZo0oX17bV+ewMBAYmNjOXHiRJm2VWSzI3F5h/HtznQu//YvLmbnOdsUl2ffvn18+eWXbNiwgd27d2M0Gi1/7GlpaTz88MP8+uuv7N+/n88//5z169czY8YMXn75ZUDTp5g8eTK//vorQUFBvPvuuzbvc+DAAcaPH8+uXbuIiooCIDU1lZdeeonVq1ezZ8+eYjL6s2fPZseOHWzfvp23336bM2fOlPkMy5YtIywsjD179pCSksKgQcUXSk+fPp0WLVqwe/duXn/9dV3vxVovFCimF2pL29N6EV9ZWpo9e/a06FpYHytXrtRlU1lY632WZVtFNjsSlx7DAGjRIABT0K9k5/6VIB8V/Vkeq1atYseOHXTq1AnQ+ukNGzakV69eNGvWjPh4bTFf27Zt6d+/P0II4uPjLfuRREZGWgR3xo4dy9tvv82UKVNK3ScqKoquXbsWO7d69WpGjBhBaGgoQLF/Vd9++22+++47QOuv//7774SEhNh8hvj4eKZMmcKTTz7JzTffTM+ePavxRirWC5UVaHtWRze0spTU+yzLtopsdiQu7zASI+vh0+R7Gga94GxTXB4pJRMmTOCVV14pdv7IkSPUqXNNJ9VgMFjSBoPBskFOyT+6sv4I/f1LD0JLKW2WX7t2LStXrmTTpk34+fnRp0+fYlqdJbnuuuvYsWMHS5Ys4amnnuLGG2/kmWeeKbN8eejRC9Wj7VlV3dDKYEvvsyzb9OqROgKX75IUcfzsVRUmXgH9+/dn4cKFnDp1CoCzZ89amuF6OHbsGJs2bQLgiy++oEePHpW691dffWXpbpw9exbQRIODg4Px8/Nj//79bN68udx6MjIy8PPzs2ziU3JDpcDAQJsfqi1K6oUWzVJY64Vaa3vm5uayYMGCYmrd5emG7t69u9RRFWchy9D7LMu2imx2JC7fwgDIu9SGnq+t4YcHupMY6XoiL7YI8w+z6HnYq76KaNOmDS+++CI33ngjBQUFmEwmZs6cSePGjXXdIzY2lrlz53LPPffQqlWrYtOQFdG2bVuefvppevfujdFopF27dsyZM4dBgwbx/vvvk5CQQOvWrUt1ZUqSnJzM448/jsFgwGQyWbZPHDx4MB9//DFhYWF0796duLg4brrpJl5//fVieUUU6YVu2LAB0Lo6RWM1cE0vdMCAARZtT7PZzJ133mnR9oTqa2mOGTOGtWvXkpmZSUREBM899xyTJk0qZnOR3md8fLxli4KXX36ZwYMHl2lbeTY7lKpOr1T3qMy0atuPO8uPfjkoT17M0n1NTeMKmp7VoSrTle5MRXqhRbi7bqhH73xWFsLrKnf1bO5sMxQehLVeaNHerSVxJS1NV8FtxjBy8s38/Ntpzl7JdbYpHkl0dHSZkYieyp133lmmswDX0tJ0FdzGYRw8dYUJs7eycu9JZ5uiUNRa3MJhhPmHMeqn7vhGzuL5lCEMXDiw4osUCoXdcYsxjJKLruw5+6BQKPTjFg6jiCs5+Xy1/TjmrJoJUlEoFMVxiy5JEULAK0v3k385xtmmKBS1ErdyGH7eXqx/oi91Gqx2tikKRa3ErRwGaFspKhQK5+BWYxigRaZ6nbmN1v+5F+/6Wtivs9WoFIragtu1MIQQJAbfyLhWfyd5QrLLbqWY1q8/+2Ji7Xak9evv1Ofp1q1blfLKo0h4pzKcP3++TJ0OWyQnJxMVFWVZk1JEVlYWvXv3tqzFaNiwIXFxccXK5Obm0qtXL8tq3qqwbNkyWrduTcuWLZk+fbrNMgcOHCimqREUFMSbb75Z7vV66nUIVY0pr+5RmbUkJSkoKCiWdgX5upIx+3tbx9i3fjvXZw88QaJPyvJl9Koj0VcZKT3raxo1aiSPHDmiJPrsRZHuglktdy+FkugrTXkSfVC+jF51JPqqIqW3atUqWrRoQVRUlJLosycvL9nH8Pcqr1TtySiJPtuUJ9FXEdWR6KuKlN6CBQsYM2ZMudc7U6JPl8MQQgwSQhwQQqQJIUrpmQshmgoh1gghdgkhfhVCVF1AQCfXNQqka/P6qpVhhbVEX1JSEqtWreLQoUMAFok+g8GgW6Jv/fr1Nu9TFYm+xMREunbtapHoK4v4+HhWrlzJk08+ybp166hbt26V3weUlugDikn0VYS1RJ81egR0ZCWl9HJzc1m0aBEjR44s9/rK1mtPKpwlEUIYgZnADUA6sE0IsUhKudeq2D+Br6SU7wkh2qBt0BztAHstjOgQ4cjq3RKpJPqKoUeiTw9VleirrJTe0qVLad++PY0aNSr3eleX6OsMpEkpD0kpc4EFwNASZSQQVPi7LlAj0xZSSg6ertxGOJ6Mkugrjh6JvoqojkRfZaX0vvjiC0t3pLzrXV2iLxw4bpVOB7qUKDMN+EkI8RDgD9gUNhRCTAYmA5hCqq8APnfjEab9uBf/ltVrtjoCU1gY+2Ji7VpfRSiJvqpJ9JUlowfVk+jz8vIqU0qvpM1Xr15lxYoVfPDBB7qud1mJPmAk8LFVehzw3xJlHgUeK/x9PbAXMJRXb3WmVYs4duaK/GLLUdl2Vodq11VdlESfe6Ek+hw3rZoORFqlIyjd5ZgEfFXogDYBPkBo1VyYfiLr+zG6c1OEMcfRt1J4GNYSfWWhJPpKo8dhbANaCSGaCSG8gdHAohJljgH9AYQQsWgO47Q9DS2L7DwzeRfjOX72asWFFWWiJPpKoyT6SlPhGIaUMl8I8SCwHDACs6WUqUKI59GaNouAx4CPhBCPoA2ATixs+jic81fzyD5xO/0/eZI6oWsBtbZEoXAUuhafSSmXoE2VWp97xur3XqC7fU3TR+O6Pix+qCexTf6C0aBN6ylFLoXCMbjdalVbxIW7xiyJLCMWQaFwBo5o5LttaHhJZq0/zHtrDzrt/j4+Ppw5c8Yh/5MUisoipeTMmTOlAs6qi0e0MAD2HD9PVl7ZI96OJiIigvT0dE6fduxYb97Jk2A9sm80YiqMDFQorPHx8SEiwr4R0R7jMP49KhGT0XkNJpPJRLNmzRx+n323DCd2/75r6ZjYYmmFwpF4jMMochbmAllsI2Q1Y6JQ2A+PcRgAP+w+wYv/28eqxxYT5KOFnqsZE4XCfnjMoCdAs1B/erYMJSvXeWMZCoUn41EtjISIerxxW5KzzVAoPBaPchhFnDifhQDC6vk62xS7kNavP3kZ2vIdPatWFQpH4XEOIyvXzA1v/Mzw9uG8OMwzxi/yMjLUTIjCJfA4h+HrbWTGyETiXST6U6HwJDzOYQAMjm/ibBMUCo/EIx0GQMqJC6zef8rZZigUHoVHTatas/nQGT74+SAF+ZXfXUuhUNjGYx3GmM5N2fr0AAxeSiRYobAXHtsl8a9z7dHUsnOFwj54bAsD4FJ2HleP3s38LcecbYpC4RF4tMMIqOOFMGbhYypbt9HtOJkKuVecbYWiluKxXRLQdu7yjfiMER2etJwbuHAgGVeuiZ671WrWc0fh/Z7Q8zHo97SzrVHUQjzaYRQhpWRP+gWSIuuRcSWD5AnJljy3Ws0aHAVtb4GNb0O7sVpaoahBPLpLUsS8TUcZNnMDv53Ut8WeS3PDc4CAFVXbb1ShqA61ooXx18Qwgny9iA4pvYmw21E3Ano8AmtfhiO2d1dXKBxFrWhh1Pf35pZ2EXh7ecjjdnsIQlpC5m/OtkRRy/D4FkaRXJ+UkH+hA/Xq9HO2SdXH2w/u3wxGE/C6s61R1CI83mFYz4CMm7XFcwK4jJoEYUBYNlzJBP/SW9la62iApqXRcvWqGjNR4XnochhCiEHAW2hbJX4spZxuo8woYBraVol7pJS329FOu/DfMe2o62sqds6tBYPPHyeix1n46V9wy3tAabGdkgrjCkV1qNBhCCGMwEzgBrSd3LcJIRYVbo9YVKYV8BTQXUp5TgjR0FEGV4d6ft4A5OSbMQiByWgo5iDcaooVoF4kZ/YFEGr4HJLGQLNeSmxH4VD0jAJ2BtKklIeklLnAAmBoiTJ3AzOllOcApJQuu678xPkser22hu92nXC2KXYhc28gBEfD4kcgL9vZ5ig8HD1dknDguFU6HehSosx1AEKIDWjdlmlSymUlKxJCTAYmA5hCTCWza4Swuj70j21E89DSU6zW3ZOitKt3UbwahXPsu0M07XOW02NaYAqLcbZJCg9Gj8OwNUpYcgNRL6AV0AeIANYJIeKklOeLXSTlh8CHAL7NfJ2yCakQgpdvsd31KOkc3KGLYhnEXPI4DYZ3okHCKOcapPBo9HRJ0oFIq3QEkGGjzA9Syjwp5WHgAJoDcVkuZufx2eajFBR4yObJg18H5SwUDkaPw9gGtBJCNBNCeAOjgUUlynwP9AUQQoSidVEO2dNQe7Nq30n++X0Ku9PPV1zYXZAStnyoHQqFA6iwSyKlzBdCPAgsRxufmC2lTBVCPA9sl1IuKsy7UQixFzADj0spzzjS8Ory14QwWjQIICGinrNNsS+H1kLaSmjeBxpc52RjFJ6GrlhpKeUSKeV1UsoWUsqXCs89U+gskBqPSinbSCnjpZQLHGm0PfAyGizOwuwp3RIh4Ob/aJGg398L5nxnW6TwMDxkcUXVmb/lKIPfWkdufoGzTbEPgY1g8Aw4sUNbBq9Q2JFa7zDC6/kS0ySQq7ke9K9x3K0QOwR+fhUun3a2NQoPwuPXklREn9YN6dPaJQNTq05R1+RkCgQ0cLY1Cg+i1juMIk6cz+K3Py/RN8ZDnId/qDbwCXDmIIS0KFVELU5TVBblMAp5blEqu4+fZ8PUfpiMHtRTO7gaPrsVRs4tlVVy3YlanKaoCOUwCvm/wbF4GYVnOQuAqB7QJBEWPYiXXx1nW6Nwczzs66g60aH+RAT7AXjOjAmAlzfcOgsKzIRffx7Mec62SOHGKIdRgmd/SGHS3G1I6SGxGaCNX9z8Jn4NcmHlNGdbo3BjlMMowXWNA0mIqOs5wVxFJIzk7AF/CGxSZhFTWBj7YmLZFxNLWr/+NWicwl1QYxgluKOL5+71cXJXXep3e1BL2GhBWc+QqAFQhS1UC6MMdh8/T96FJGeb4Rh+XwGzByK8PGisRlEjKIdRBu+tTSMnsz/5Zg/8qIzekL6NsC7nocADn0/hMJTDKIMXhsXhH/0OXp42zQrQvDfc8AJBkdmw9hVnW6NwI9QYRhk0DPQhPCiEuDnxyLz6GLzPuoVkn26uf4Dz771EPV6D0OsgYWSx7KIBUOu0igJVeOA/n/Zj+YjljKi/AN+T/2T9qJ3Fdn13e4Tgjx11tcCuYxtLZbdcvYrY/fssh3UIuaL2oloYFTCqYyQxjQMJrOOBr6pAwB1fg8nX2ZYo3ATVwqiANmFBjO7cFINB2JqJdH+8/bTVrZlpMH8kXD3rbIsULoxyGDrZkJbJ1aP3cSnbQ0OrL/+pyft9fhvkXnW2NQoXRTkMnfiYDFBg4uyVXGeb4hiie8DwjyB9G3wzScn7KWzigR1zx9Ahqj4t4xZy8+JrsnceNWsC0HYYXH4Nlj4O/3sU/vqW1l1RKApRDqMS/DRyOTn5Zt5dc5CxXaPo921HZ5tkf7pM1ronR9ZDXpY2xqFQFKIcRiVJP5fF+z8fJDTQg7Ul+v0L8rO12RNzHhids62lwvVQDqOStGgQwKrHehMR7MdrB5xtjYMQQnMWuVfh81HQop+zLVK4CMphVIEioZ2C3BD2ZlykTViQky1yEF51IKARrHqO4Ouq/oxKO9RzUNXNUbIAABbxSURBVA6jihQUSLLSxzH121/54YHuDPpmULFIUI8YEDUY4Zb3wZxDY36Ejf+Fbg9VuhqlHeo5KIdRRQwGgU+Tr3j3tu8QQpBxJYPkCcmWfHfY+V0XRhOM+ISLd0cQ9NM/weAFXe9ztlUKJ6ErDkMIMUgIcUAIkSaEmFpOuRFCCCmE8MDpg9IYfTMs3RNzjgfv/2E0cWJTMHS5D1ooJa7aTIUtDCGEEZgJ3ACkA9uEEIuklHtLlAsE/g5scYShrsxnm49y9dA/SDlxgbjwuoDWJbFuZbh7F8XUJJx9j/wA/ABI6iX40mTBDjCo2L/ahJ4uSWcgTUp5CEAIsQAYCuwtUe4F4DVgil0tdAOGJIXx/MbXiG3yV8u5ks7B3bsoxQYpf18B80fAd5Nh6LuaMrmiVqDnn4dw4LhVOr3wnAUhRDsgUkq5uLyKhBCThRDbhRDbzZfMlTbWVQnyMVEn5BeMBsGFrDwu53h4WHXLAZzaEwjJX8MXt0HOZWdbpKgh9LQwbMUGW9ZtCiEMwH+AiRVVJKX8EPgQwLeZr8et/cwzFzDq/U00b+DPe2M7ONscxyEEFy+0Jn/r7zQpWE32lGiO/1Ifc44RUNOmnoweh5EORFqlIwBrNZVAIA5YK7R1B42BRUKIIVLK7fYy1BWxHqcI8w/DZDRwZ49ookP8nWyZ47E4hANL8f3mbq779gNo0RdQ06aejB6HsQ1oJYRoBpwARgO3F2VKKS8AoUVpIcRaYIqnOwsoPU4BcFunppbff1zIokld1xGnsSW7V21a3wT/+BX86mvpC+nVr1PhslToMKSU+UKIB4HlgBGYLaVMFUI8D2yXUi5ytJHuyIa0TP72yTY+GN+Bvq1dY0d4h3UTipzFwTUwfyTBLf20fU90rHRVUaDuha7ALSnlEmBJiXPPlFG2T/XNcn86RAUzsXs0HaOCnW1KzRHRCVoOoHHBUlj0EAyeASafci9RUaDuhZpEdxA+JiP/NziWQB8T+eYCCnJrgeOoEwCj55OZGgC7PoU5g1UXxcNQDqMG+NcPqVw98gDnr3qoWpc1BiOnk4Pgtvlw+jf4bZmzLVLYEeUwaoDJvZrj3WA59fxqUYBT7M3w4FboOElLn9wLBZ4Te1NbUYvPaoBmof54B28DIO3UJbwMBqJDPX/qlSBtFsboY4ZZN0J4Oxj+scNupwZQHY9yGDVIQYHkvs92EujjxTf3dUPUEr1Mc7YBbpoO/5sC73fHv7H+lkZlnIAaQHU8ymHUIAaD4K3R7fD2Eh7tLErHe4RDu7EQ3gG+/htN++yDJY/DoFcrrEs5AddCOYwaxlqda9b6w3SICiYpsp4TLbI/ZXYDGsbC5LWcGdeckPbZNle6OiS4TGE3lMOoIUqGkX/71/8xd+MR0k5d8jiHUS4mH07trkvIF9p2DT7BubDqBej9BHjVUWMOLo5yGDWEdRh5/Nx4/Ot48e393ajrqyly55kLMBlr0aRVYZcsICwH1s2A/YthyDsQ2Ul3FbbGNxSOpRb9hboeoQF1MBkNXM3NZ8R7G/l43SFnm1TjZKYGwh0LtSXys26AZU/pXi5fNL5RdKjWieNRDsMFMAhBdKh/rVjlapNWN8D9m6Dj32Dzu7BjjrMtUpSB6pK4AD4mI2+NbmdJbzyYSZsmQR4b6GU9sGnpRvgEwc3/gaSx0DhOO3d0EwQ2hvrNnGSpoiTKYbgYl7LzuPfTHfSLacibVk7Ekyi36xBRKDxUUAA/PgznDkOXe6DnFPCtRYPDLopyGC5GoI+JetGfs+LKb8TPvYyUgq+t8gcuHOh5+5/YwmCA8T/A6hdh4zuwaz70eQqEfYTaVFRo1VAOwwlUpCh+TuwkdVIyUkoeXrCb2W0+4/XCPI/d/8QWQU1g2EythbH8/2Dp4wQ0qW+XqlVAWNVQDsMJlGwRDFw4sJQDAcgvkAT5emHMz65R+1yOJgkw4Uc4sp7LC+7WzqV8Q516eXa7ha2AMdXiKI1yGC5AWV0Kk9HAi8Pi2Tt1NfFz4zFnRVDP2LOGrXMRhIBmPQGh7Si/chrNB52Gr8ZD7yehUdtqVV/SOagWh23UtKobIIDkCcn0DXwJcWYE2Xm1fJm40QT3/MLplABIWw3vdYMvx0JmWrFiRa2GfTGxpPVTO7bZA9XCcCP+c1sSJ85n4WMyUlAgScm44GyTnIdvMJkpQTSYtwE2vwdbPoCel7S8AjMYjMVaDarFYB9UC8ON8PYy0KxQR+OLbccYOnMD5qwIJ1vlZPzqQ7+n4bF9EFY4Db3o77DgDjixw7m2eSCqheGmDEsKR0qYvl/TzMy8nENoQB0nW+VEvK2iZOtHw4b/autTIrvC9Q9gCmuiVsHaAdXCcFP863gxtmsUQsDZK7nc+J9feGf17842yzXo9Tg8kgKDpsOlP+CrcbT81w1q3YkdUC0MD8C/jpGxXaMY0KYRAxcO5MSlUyAKEKKgWDmPDfKyhU8QdL0POk+G/f+DsCTt/JENcGCJdj44yrk2uiHKYbg5Yf5hdJyvfQyfZGjp0Q3ms+HgGb67vxs+JqOlrCcEedlch1IeBiO0GXItnbFTGyTdNBNa9of2E7Td24wmB1nsWSiH4ebYajEsT/0TvzpeFmdx/mquxyxkq3ZXottD0PYW2DEXdn0GX43TBksnr7WHeR6PLochhBgEvIW2VeLHUsrpJfIfBe4C8oHTwJ1SyqN2tlWhk4FtGzOwbWMAjmReYdBbv/DqrQlOtsqFqBuhzaz0fhLSVkJuof6GOQ++mQSxQ/COaFysJaPGPDQqdBhCCCMwE7gBbSf3bUKIRVLKvVbFdgEdpZRXhRD3Aa8BtznCYEXlCPTxYnSnplzfPAT2aBtE+5qMHtPiqBZGL2g96Fr6/DHI2AV7f6BFv0CtKxM/gn03318j5rjDgjg9LYzOQJqU8hCAEGIBMBSwOAwp5Rqr8puBsfY0UlF1QgLqMG3ItbDpFxbvZdex8/zyRF9u/u6m2rHyVS8hLeDve+Dwz5C8EPYtgt3zCYxpzb6YWIRRIs2aCrojPmR3WBCnx2GEA8et0ulAl3LKTwKW2soQQkwGJgOYQtQgkzP4e/9W/H7yMiajgYwrGTxx3Y8MaNOIJnV9iy2Cq7XOw2CAFn214y8zIG0VEc8M1s4vnQq//0Tmuv3wxx5onKBrh3pPQo/DsPVGbIoSCCHGAh2B3rbypZQfAh8C+DbztY+wgaJSxDQOIqaxttVBQW59pv24l5z8Au7q2byUUHGtx+SrbflYRGRnOJlCSOxB+KAX1G0K7e6APlOdZ2MNo8dhpAORVukIIKNkISHEAOBpoLeUMsc+5ikcicH7LGun9CHYXxvPWLP/FJ9uPsr0W5WzsEnccIgbzu+J13HdJ/+EfYvh0p9anpSwchpEdddW1Zp8AfcYl6gMehzGNqCVEKIZcAIYDdxuXUAI0Q74ABgkpTxldysVDiOyvp/l94WsPE5dyibYz5sw/zDafNgPg+k8wpBXe7soNjDnGKH9eO2QkrR+/ZHnjtN88CmMpjcpyIerp+pw+Y86GL2jaeni4xKVoUKHIaXMF0I8CCxHm1adLaVMFUI8D2yXUi4CXgcCgK8LtwA8JqUcUmalCpdkWLtwhiaFIYRg2a3LGPjmL4QG1OHzu7uqLooVtsR2Wu36DfJz4Mg6DL+vJCBtBQFhaTDyMa3Q+eNwej/CWDz61roF4g7rW3TFYUgplwBLSpx7xur3ADvbpXAS1nu+vjA0DrPUhppkgZEJs7dyT+/mdGsR6izzXIIyuxRedaDlAO1gOpw9BAGNtLyUb2Dls1w3HPhkMDTrBc16kffHCWL3768p06uNivRU2EQIQZfmIZa0zAsm/dJV8s2aAzl9KYffT12iS7MQjIbyZwpqjXBxSeo3v/a7yz3QOI4Lr9yPz4Wt+BzZAKtewTu8cEn+0Y1g9IYmiS4dpq4chkIXhjqZrLz72uTX97tO8NKSfayd0ofoUH+ycs34mAw2d6WvVcLFZWHyhZYDCJ71m5bOOgcnU2nxQg8tvep5OLYJvHwhoiM07YpfQ9ebO1AOQ6GLMP8wEuZdCy9v7NOUOX/7hOhCQZ8X/6cFhC1+qAeGClocCsA3GKJ7XEuPmgdHN8CxLZrjWPcGITFWn+faV7XVtWHtIKSVFhdiRUWzMSXzq4pyGApd2FI6f2hzfy2uF8i7mEBBfn0SP71f20ZBvkZskyAm9VC7lukioKG2KK7tLVo65zKZwwZxPCYWYSyg1dCTGL0LQ5e8A7Xl+p0maeWlJC+j+FhIydmYYlGk1Qg2Uw6jFmFrP5SqUt4YRNycBMLJ5VK2tg2AlPDsDykMSQqnQ1RwsbK1dnyjIuoEEL10/bV0gZmDXVvT4oPntCX6J3ZC9kUt7+whWt1yEj4dDuHtIawdXr5m7cXbORJVOYxaRE19iEJIZk/oZEnL/Hp8vzuD+Ih6dIgKRub78d7agwxrF6bGN/RiMCIDotg35kWrk9OB6ZgC8mnQoT51L5+CdW+ANNNqKJpQUMxfIDONei2uwPFt0LB6cSDKYSgcjsF0nh3/HGCZojVnR/Lqsv10aa7tYnbw9GW2HT7LzYmlWzyVaYF4emtFV4Ro7lX4M5k/7x9J48c6aucOraFJpwswq/rRD8phKOyOra6Pl9Fg+WPzCjjA1qf7U79wif3KvSeZvmw/N7ZtTJh/GLHvDkHmheBVdyfhAY11t0BUawXw9oOmXTiX5k/jwMIYkI6TSLvvNVrOfwtOpsJzT1S5euUwFHZHz7/qDQN9LL8n92rOwLaNqe/vzfIRy/m/75L5KfVPtv19NkIIvt2ZztVcM2O7FtfgtNWiUNjAYCDvihfEDNYOlMNQuBHWLZAwfy0UvWh6FuDFoXH8vV8rS0zH0pQ/OX811+IwXlmyj9CAOqVaFArHoxyGosapqAViMAga173WAvlofEeu5uZb0gdPXyHLarvIYTM30C+mIX/v3wqAE+ezaBJ07XpbePp4h631LvZAOQyFW+Dnfe1P9eMJ2mDed3PBXCBp3SjQ4mBy8s30fm0N9/dpQZh/GHFzEsg735HwBhdZfcdCSx0lWyfW4kHg/g7EUUvolcNQuDVGg+DVEdciUAsK4KVb4mgbVpdHw5dz8PRl+v/7Z676rCB+bjwFeUFk/3Er4ZHtAc3BXM0xl3IOtXLAVAfKYSg8Cl9vI7d1ampJNwvx55fH+1LX90bq+pnY98dFHl+4h+f6TgRgx9Fz3P7RFj6/uwvdWoRyJPMKaw+coiDfr4w71G6Uw1C4FSUHTCvCYBA0Dbn28cc2CWLxQz0t6chgP54eHGuRLdx65CzTftxLs7go4ufGk3chCXPmTax7bDiN6/rw28lLHDx1mX6xDanjZSx1P09HSOkcaU3fZr4y63CWU+7tbuyLiS2mJq1wHFJKTl/KISSgDkaDYGNaJvd8OxfZ8FOEMJNzuj+5mQPY/8JN+JiMfLLhMN/vzuDb+7phNAh2HTvHyYs5DIprXOl719RArBBih5SyY1WuVS0MhcIKIQQNrWZYurUMJfmJxwBNOetyTj6d5wzCx6SJAwfU8SKinq9FE+Sr7cdZsfeUxWEkvfESFy8F4x/9HgCBuT2Y0uFphrePADRZRD9vo0XF3dUDz5TDUCgqQUAdLyJDjKU+5vi52n+l2YeChoHEz30EAD/fAdzT5laeHKRthtTihX/zyYYjFodx/X8+IjvXgH+zmYT5hzFj+QF8TAYe7KdNEW86eIZAHy/iwutq9UtpU3OkplAOQ6GoJNXpJrSIWcHhSwuIn5sNQL0GPXmy09MMTboXgIe+2IWv6ZrWxbRFqUSH+vHBOK0HMejNdSRE1OX1kYkAPPdjKrFNgph15C4yrmSQf/k6hOkckfVNLB+xnNz8Ary9DKW6O1VFOQyFogb5aeSycvP/O6ad5XeYfxjHc6aRniWJn5tJmH8Yw9u/Wyyobdex8/iajJbuTPyzyxnRMYKFZ0cDkPjcT4zvFkXGlQx+Hf8rY2dtIYXrq2y/chgKhYtiS7TonSPDAPjnnsKTAXDw5LUZo6/uvZ6AOl4s/FHrvjzYryWJEfX4/BTkmSV55upNciiHoVC4CXq6QrFNgiy/hRA80LelltgI3l4GvrrnesS9VbdBOQyFwgOxp7qaNcphKBQeiKPWwRgqLqJQKBQaymEoFArd6HIYQohBQogDQog0IUSpve2FEHWEEF8W5m8RQkTb21CFQuF8KnQYQggjMBO4CWgDjBFCtClRbBJwTkrZEvgP8Kq9DVUoFM5Hz6BnZyBNSnkIQAixABgK7LUqMxSYVvh7IfCOEELIcla2NcvLh5fDi59s2hXGfqP9fr8HnD1cPL/lABhVGIP7dju4fKp4fpthMGym9nvGdZB7pXh+0u0w+HXtd8l7A3SeDAOehZxL8O+Y0vk9HoFeU+DSn/DfDqXz+/0Lut4LZw7CB71K59/0KrQbC3/s0TbkLcmQ/0Lc8NLnFQoXocLVqkKIEcAgKeVdhelxQBcp5YNWZVIKy6QXpg8WlsksUddkYHJhMg5IsdeD2IlQILPCUjWPK9qlbNKHK9rUWkoZWJUL9bQwbK10Kell9JRBSvkh8CGAEGJ7VZfYOgpXtAlc0y5lkz5c1aaqXqtn0DMdiLRKRwAlV7FYygghvIC6wNmqGqVQKFwTPQ5jG9BKCNFMCOENjAYWlSizCJhQ+HsEsLq88QuFQuGeVNglkVLmCyEeBJYDRmC2lDJVCPE8sF1KuQiYBXwqhEhDa1mM1nHvD6tht6NwRZvANe1SNunDo2xymkSfQqFwP1Skp0Kh0I1yGAqFQjcOdxiuGFauw6aJQojTQojdhcddNWDTbCHEqcKYFlv5QgjxdqHNvwoh2ruATX2EEBes3tMzNWBTpBBijRBinxAiVQjxsI0yNfqudNpUo+9KCOEjhNgqhNhTaNNzNspU/tuTUjrsQBskPQg0B7yBPUCbEmXuB94v/D0a+NIFbJoIvONIO2zY1QtoD6SUkT8YWIoW89IV2OICNvUBFtfwe2oCtC/8HQj8ZuP/X42+K5021ei7Knz2gMLfJmAL0LVEmUp/e45uYVjCyqWUuUBRWLk1Q4HCeG8WAv2FY2WR9dhU40gpf6H82JWhwDypsRmoJ4Ro4mSbahwp5R9Syp2Fvy8B+4CScf41+q502lSjFD775cKkqfAoOcNR6W/P0Q4jHDhulU6n9Iu0lJFS5gMXgBAn2wRwa2FzdqEQItJGfk2j1+6a5vrCZu9SIUTbmrxxYRO6Hdq/ntY47V2VYxPU8LsSQhiFELuBU8AKKWWZ70nvt+doh2G3sHI7oud+PwLRUsoEYCXXvLAzqen3pIedQJSUMhH4L/B9Td1YCBEAfAP8Q0p5sWS2jUsc/q4qsKnG35WU0iylTEKLzu4shIgrabKty8qr09EOwxXDyiu0SUp5RkqZU5j8CLCxNLXG0fMuaxQp5cWiZq+UcglgEkKEOvq+QggT2oc5X0r5rY0iNf6uKrLJWe+q8H7ngbXAoBJZlf72HO0wXDGsvEKbSvR3h6D1SZ3NImB84QxAV+CClPIPZxokhGhc1OcVQnRG+3s64+B7CrTI4n1SyjfKKFaj70qPTTX9roQQDYQQ9Qp/+wIDgP0lilX+26uB0drBaKPGB4GnC889Dwwp/O0DfA2kAVuB5i5g0ytAKtoMyhogpgZs+gL4A8hD8/yTgHuBe+W1Ue+ZhTYnAx1dwKYHrd7TZqBbDdjUA63Z/Cuwu/AY7Mx3pdOmGn1XQAKwq9CmFOAZG3/nlf72VGi4QqHQjYr0VCgUulEOQ6FQ6EY5DIVCoRvlMBQKhW6Uw1AoFLpRDkOhUOhGOQyFQqGb/weDsrNNkGqA5AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.plot(\n",
    "    bins,\n",
    "    lamb(bins) * np.exp([-sp.integrate.quad(lamb, 0, x)[0] for x in bins]),\n",
    "    ls=\":\",\n",
    "    label=f\"$p(t;\\lambda(t)={background}+{base}^{{-t}})$\",\n",
    ")\n",
    "\n",
    "pl.plot(\n",
    "    bins,\n",
    "    np.where(\n",
    "        bins < cond,\n",
    "        0,\n",
    "        lamb(bins)\n",
    "        * np.exp(\n",
    "            [\n",
    "                -sp.integrate.quad(lamb, cond, x)[0] if x > cond else np.nan\n",
    "                for x in bins\n",
    "            ]\n",
    "        ),\n",
    "    ),\n",
    "    ls=\"--\",\n",
    "    label=f\"$p(t|t>{cond};\\lambda(t)={background}+{base}^{{-t}})$\",\n",
    ")\n",
    "\n",
    "pl.hist(\n",
    "    x,\n",
    "    bins,\n",
    "    density=True,\n",
    "    fill=False,\n",
    "    histtype=\"step\",\n",
    "    label=f\"empirical s.t. $\\lambda(0)={lamb(0):.3f}$\",\n",
    ")\n",
    "\n",
    "pl.hist(\n",
    "    x[x > cond],\n",
    "    bins,\n",
    "    density=True,\n",
    "    fill=False,\n",
    "    histtype=\"step\",\n",
    "    label=f\"empirical s.t. $\\lambda({cond})={lamb(cond):.3f}$\",\n",
    ")\n",
    "\n",
    "pl.legend(loc=\"upper right\")\n",
    "pl.xlim(xmin=0, xmax=xmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fit GMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loglik=-2.531 elapsed=0.0s\n",
      "2 loglik=-2.225 elapsed=0.2s\n",
      "4 loglik=-1.799 elapsed=0.5s\n",
      "8 loglik=-1.577 elapsed=0.9s\n",
      "16 loglik=-1.495 elapsed=1.8s\n",
      "32 loglik=-1.467 elapsed=3.5s\n",
      "64 loglik=-1.458 elapsed=6.8s\n",
      "100 loglik=-1.452 elapsed=10.4s\n"
     ]
    }
   ],
   "source": [
    "model = GMMModel(x[:, None], num_clusters=num_clusters)\n",
    "trainer = GMMTrainer(model)\n",
    "\n",
    "for t, epoch in elapsed(range(100)):\n",
    "    trainer(x[:, None])\n",
    "    if (\n",
    "        np.allclose(np.log2(epoch + 1), np.round(np.log2(epoch + 1)))\n",
    "        or epoch + 1 == 100\n",
    "    ):\n",
    "        loglik = model(mx.nd.array(x[:, None]))[0].mean().asscalar()\n",
    "        print(f\"{epoch+1} loglik={loglik:.3f} elapsed={t:.1f}s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inferred lamb(0)=1.366, lamb(1)=0.908\n"
     ]
    }
   ],
   "source": [
    "lamb0 = infer_lambda(model, xmin=0, xmax=1)\n",
    "lamb1 = infer_lambda(model, xmin=cond, xmax=cond + 1)\n",
    "print(f\"inferred lamb(0)={lamb0:.3f}, lamb({cond})={lamb1:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mixture_pdf(bins):\n",
    "    log_marg = model(mx.nd.array(bins, dtype=\"float32\"))[0]\n",
    "    return log_marg.exp().asnumpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1.3621529662981628)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARsAAAD5CAYAAAAX1pEcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deXhURda438pGQhIIO0IwAdlCSIAIyCogqAh84IKOCIOyiDI6P8dlFEdHHPUTZnEcxWVkXBAX8AOVYUQZAXFEAWWXfU2AAA4QFlmy5/z+uN2drdPdgXTfTjjv89yn771Vt+rk0n2oOnXqHCMiKIqi+JsQuwVQFOXSQJWNoigBQZWNoigBQZWNoigBQZWNoigBQZWNoigBIcxbBWPM28Aw4KiIdHRTPhp4zHF5FpgsIpu8tRsSEiJRUVGVFFdRlIvl/PnzIiIBH2gYb342xpirsZTI7AqUTS9gu4icNMbcADwtIld56zg6OlrOnTt3gWIrinKhGGPOi0h0oPv1OrIRkW+MMYkeyleWuFwNxF+8WIqi1DSqeig1AfiiittUFKUG4HVk4yvGmAFYyqaPhzqTgEkAERERVdW1oijVgCpRNsaYVOBN4AYRyaqonojMBGaCZbOpir6V4CQ/P5/MzExycnLsFuWSJTIykvj4eMLDw+0WBagCZWOMuRz4BPiliOy6eJGUmkBmZiaxsbEkJiZijLFbnEsOESErK4vMzExatmxptziAb0vfc4D+QENjTCYwFQgHEJG/A08BDYDXHF+qAhHp6i+BlepBTk6OKhobMcbQoEEDjh07ZrcoLnxZjRrlpXwiMLHKJFJqDKpo7CXY3n/weBBrXB3FBhYuXMj06dM91pk1axaHDx8OkETu+frrrxk2bBgAubm5DBo0iM6dO/PRRx/ZKldlqLLVqItiWgtokgzjF9stiXKJMXz4cIYPH+6xzqxZs+jYsSPNmjXzud2CggLCwvzz89qwYQP5+fls3LjRL+37C9tGNjn5hSROWUTilEWkZ0eRfjDTLlGUGkhGRgbt27dn4sSJdOzYkdGjR7N06VJ69+5NmzZt+OGHHwBLkdx///0AjBgxgtmzZwPwxhtvMHr0aObPn8/atWsZPXo0nTt3Jjs7m8TERI4fPw7A2rVr6d+/PwBPP/00kyZN4rrrrmPs2LEUFhby29/+lm7dupGamsobb7xRoZx33nknqampjBw5kvPnzwOwePFi2rdvT58+ffjkk08AOHr0KGPGjGHjxo107tyZvXv3+vU9VikiYsthwmqJi7mjZdfvk0SpOWzbts3W/tPT0yU0NFR+/PFHKSwslLS0NBk3bpwUFRXJggULZMSIESIi8s4778h9990nIiI//fSTXHHFFfLNN99ImzZtJCsrS0RE+vXrJ2vWrHG1nZCQIMeOHRMRkTVr1ki/fv1ERGTq1KmSlpYmjr1H8sYbb8izzz4rIiI5OTly5ZVXyr59+8rJCci3334rIiLjxo2TP//5z5KdnS3x8fGya9cuKSoqkltvvVWGDh0qIiLLly93nXvD3b8DcE5s+M0Hh81m+79oE3LIbimUGkbLli1JSUkhJCSE5ORkBg4ciDGGlJQUMjIyytVv0qQJzzzzDAMGDOCFF16gfv36le5z+PDhODcYf/nll8yePZvOnTtz1VVXkZWVxe7du8s906JFC3r37g3AmDFj+Pbbb9mxYwctW7akTZs2GGMYM2ZMpWUJNoLDZuOkMB9Cg8MBSan+1KpVy3UeEhLiug4JCaGgoMDtM5s3b6ZBgwYeDcJhYWEUFRUBlHNajI4u3t8oIsyYMYPrr7/eo5xlV42c18G2mnSxBMfIxkn+ebslUC5hfvjhB7744gs2bNjAX/7yF9LT0wGIjY3lzJkzrnqJiYmsW7cOgI8//rjC9q6//npef/118vPzAdi1axfuIh0cOHCAVatWATBnzhz69OlD+/btSU9Pd9lk5syZUzV/pI0Eh7IxobxSMAIi69otiXKJkpuby913383bb79Ns2bNeOGFFxg/fjwiwl133cW9997rMhBPnTqVBx54gL59+xIaGlphmxMnTqRDhw6kpaXRsWNH7rnnHrcjqqSkJN59911SU1M5ceIEkydPJjIykpkzZzJ06FD69OlDQkKCP//8gOA1no2/CAmPlKL8HCgqhGfq89f8kTz0v2/ZIotS9Wzfvp2kpCS7xQh6MjIyGDZsGFu2bPFL++7+HYI2no3fKcwD4KHw+fz6d035V1EvAJrHRfHdlGvslExRlCrE/mlUWCQ8uBWAGUObkDF9KBnTh3LoVLbNgimK/0lMTPTbqCbYsF/ZGAN1mltK59xRu6VRFMVP2D+Nyj4JK16Aghw4Gzw7VBVFqVrsH9mcPwErZ1jn51TZKEowYIx52xhz1Bjjdo5njBltjPnRcaw0xnTy1qb9ysZhIKZRErTwmpRBUSrF5s2badq06SVjF6lCZgGDPZSnA/1EJBV4FkcETk/YP41yKptrnoSkYfbKoviPd4aWv5d8I3S/G/LOwwe3li/vfAd0GQ3nsuD/xpYuG7fIp26ff/55Vq5cyZNPPsmHH354AYJfmogfsqrYr2wKHMomrFZxTJsa5qat2IfT81YVjV/xKauK/crGObLZNBfmjIKHd0J0A3tlUqoeTyORiNqey6Mb+DySUXwizBiztsT1TLGSEVQaX7KqOLHfZpPYG546AW0HQ1G+GomVKmXz5s0kJCTw+uuvl7qfnZ1Nv379KCwsBKzYMe3ataN169auyH15eXlcffXVFW7a9IXx48fTuHFjOnYsl0zWRU5ODt27d6dTp04kJyczdepUV9mpU6cYOXIk7du3JykpybWHyluZFwpEpGuJ40IVjTOrygjxkFXFhR1xLaRsPBsRkb1fi0ytI7LvPyIikvDYZx6idCjBjt3xbEqycuVK6dGjR6l7r7zyivztb38TEZGCggJp1aqV7N27V3JzcyU1NVW2bt0qIiJPP/20vP/++xfc93/+8x9Zt26dJCcnV1inqKhIzpw5IyIieXl50r17d1m1apWIiIwdO1b+8Y9/iIhIbm6unDx50vWcpzInFxPPBkgEtlRQdjmwB+jlS1sSFPFsDm+Azx4Csf6H0ZGNUtU0btyYrVu3lrr3wQcfMGLECMDa7d26dWtatWpFREQEt99+O//85z8BuPHGG/nggw8uuO+rr77aa1wcYwwxMTGAlW8rPz8fYww///wz33zzDRMmTACsxI5xcXEAHsuqAkdWlVVAO2NMpjFmgjHmXmPMvY4qJbOqbCwzLXOL/comay+sfQvCrIBD6tinVDVTpkwhNzeX/fv3A9b0aN++fSQmJgJw6NAhWrRo4aofHx/PoUNWMLeOHTuyZs2acm327duXzp07lzuWLl16QTIWFhbSuXNnGjduzLXXXstVV13Fvn37aNSoEePGjaNLly5MnDjRFaLCU1lVICKjROQyEQkXkXgReUtE/i5W+iZEZKKI1BORzo7Da/om+5VNQa71GdsEuo6Hxu3tlUepUSxevJhz584xdOhQ1+jm+PHjpUYB4ibygTNwVWhoKBEREaXi2QCsWLGCjRs3ljsGDRp0QXKGhoayceNGMjMz+eGHH9iyZQsFBQWsX7+eyZMns2HDBqKjo132JE9lwYr9yqbQoWzComDYi9Cqv53SKDWInJwcHn30UV577TVSUlJcjn1RUVGlIuzFx8dz8OBB13VmZmapTAq5ublERkaWaruqRzZO4uLi6N+/P4sXLyY+Pp74+Hiuuspydh05ciTr1693yVxRWbASBEvfVhQzQiOs2DZ55yCyjr0yKTWC5557jrFjx5KYmEhKSgoLFy4EoF69ehQWFpKTk0NkZCTdunVj9+7dpKen07x5c+bOnevyy8nKyqJRo0bl8mWvWLGiyuQ8duwY4eHhxMXFkZ2dzdKlS3nsscdo2rQpLVq0YOfOnbRr145ly5bRoUMHAI9lwYr9IxuwdnyHRVhepO/dZLc0Sg1g586dLFmyhN/85jcApUY2ANdddx3ffvstYMUUfuWVV7j++utJSkritttuIzk5GYDly5czZMiQC5Zj1KhR9OzZk507dxIfH89bbxUHiBsyZAiHDx/myJEjDBgwgNTUVLp168a1117rSkg3Y8YMRo8eTWpqKhs3buR3v/ud63lPZUGJr8tWVX2UW/oWEfn4bpEXO4qILn1Xd4Jp6dsd69evlzFjxnitd9NNN8mOHTsCIJF/qFapXHzY/WmMMS8bY/Y4doCmXbDmi25krUbZFKpUuXTo0qULAwYMcDn1uSMvL48bb7yRdu3aBVCymosv06hZeN79eQPQxnFMAl73ULc8P86DBfdZ5zGNoSDbstsoip8ZP368x4DlERERjB07tsJypXJ4VTYi8g1wwkOVEcBsxwhtNRBnjLnMZwkOb4BtlgMV0Y2sT43Ypyg1jqowEDcHDpa4znTc843C3OLEdM3SYMATEBFTBWIpihJMVMXSt7t4EG6NLsaYSVhTLQhxdF2YZ4WXAMuhz+HU1zwuisQpi1znmmlBUao3VaFsMoEWJa7jAbe5S8XaXToTrLxRgBXPJjTCWQHOHIGwyFLKxal0lOqFiNS4FLLVCQmyhZaqmEYtBMY6VqV6AKdF5IjPT0dEQ0wT67wwH/6aBD/8owrEUuwkMjKSrKysoPvCXyqICFlZWeU8n+3E68jGsfuzP9DQGJMJTAXCAcTalPU5MARru/l5YFylJBj21xLSREBknBqIawDx8fFkZmZy7JhurLWLyMhI4uO9RusMGF6VjYiM8lIuwH1VJlFMYw0zUQMIDw+nZcuWdouhBBH2b1dY9gws/UPxdXRjDTOhKDUQ+zdi7l9ZvDIFENMIftpsnzyKovgF+5VNQS7Uji6+Thtrpe5QFKVGYb+yKcwvXvoGuEL9aRSlJmK/zaYwt7SyyfkZMtdZicsURakx2K9s6jSDOiV2N2SsgDevgeM77ZNJUZQqx/5p1Nh/lr6Obmx96oqUotQo7B/ZlCW6ofWpjn2KUqOwX9m8fwusebP4urYjx072SXvkURTFL9ivbNJXwKkDxde16oAJVWWjKDbijwid9iobEcdqVK3ie8bAyLeg40j75FIUZRZVHKHTXgOxM41LWETp+8maYUFR7EREvjHGJHqo4orQCaw2xsQZYy7zFPHB3pFNYZ71GVpG2RzZBAe+D7w8iqL4SqUjdNo7spEiuKwzxJYJWbzsWTh/HCZ9bYdUilLTCTPGrC1xPdMR2K4y+Byh09VpJTuoWiLrwD3/KX8/qh5k7Q68PIpyaVAgIl0vsg2fI3Q6sX81yh1R9XQ1SlGCm0pH6LR3ZHNyP3w0BgY9Da0HFt+Pqgc5p63c3yGhpYKfgwZAVxR/448InfYqm9wz8NOPkHe29P2oetZnzmmoXb+cYtEA6IriX/wRodPmpe8KVqPaD4UmHaxg6Iqi1AiCU9nEtbAORVFqDPYaiAtyrc+wWqXvZ5+CLR/DqYPln1EUpVpir7KpFQMJfYptNE7OHoX54+HAanvkUhSlyrF3GtX8ShjnxtjrMhCfcv+Yrk4pSrXD/uBZ7oiKsz4r8LXR1SlFqX7YO43a/hm8nFY6xARAaDhExKhjn6LUIOxVNtkn4MRe3G6zUC9iRalR+KRsjDGDjTE7HYFyprgpv9wYs9wYs8ERSGeIT71XtPQNMGoOXPOkT80oihL8eLXZGGNCgVeBa7E2X60xxiwUkW0lqj0J/J+IvG6M6YDlypzotfcCh7IpG88GoGmK18cVRak++DKy6Q7sEZF9IpIHzMUKnFMSAeo4zuviZfenC08jm/QVsOED721s+ogrzCGfulMUxT58UTa+BMl5Ghjj2LD1OfBrn3qvlwBtB5cOC+pky3xY+rTn57P2wqeTeCv8Lz51pyiKffiibHwJkjMKmCUi8Vg7Qd8zxpRr2xgzyRiz1hizVooKrfCfd3wEoW5mc5Fxlp+NeIjHs+4dAOLNMc/1FEWxHV+UjS9BciYA/wcgIquASKBh2YZEZKaIdBWRriYk1HOvUfWsaVa+hzS83SZCfHfCTBGc2u/Dn6Ioil34omzWAG2MMS2NMRHA7ViBc0pyABgIYIxJwlI23lNaLnsGZlQQMMzpRexp+bteIox8i8fy74bIul67UxTFPrwqGxEpAO4H/g1sx1p12mqMecYYM9xR7WHgbmPMJmAOcJcj3oVnsk9WrExcXsTutyyw9GnLiBx3OR8VDii/v0pRlKDCp+0KIvI5luG35L2nSpxvA3pXuveCvPI7vp1cMRB+s7l8MHSAo9vh2xctBdOyL53rnuPXv3uCfxX1AnSvlKIEIzancsl1v+wN1o7wuMutrQtlWfuO9Vzn0QAs6H+MGRGvkPG7LmRMH8qhU9l+FFpRlAvB/rxRFSmbvHOw4gXIXFvm/nnYNBeShkO0wwadYI1o2L/Sf7IqinJR2KtsWvSAdjdUXL7sGcj4tvS9rZ9A7mnoOr74XpMUiIiF/d/5R05FUS4ae0NM9PxVxWXhtSEkvLwBOT8bOtxYPJoBy0/n8h46slGUICY480YBGONI6VJmNar73XDbu1Z5SRJ6wbEdcP5E4GRUFMVn7B3ZzBoGtWKtHd7uKBtmoqjI+gxxoyPT7rQMxrXrl4rkpytTihIc2Duyyf3ZSkRXEWWVzX83w/PNYM+y8nWjG0BsE8CK5JcxfaiuTCnKBeKPsDI2p3LJd7+07eSOjyzbjZNjO6EgG+o0c19/2z/hwPcw+PmqlVNRLiH8FVbG/lQuFTn1geVFXDLWzbGdYEKh/hXu6x/dDqtfszJpKopyofglrEwQ+Nl4UDa7l8K/nyi+PrYD6rdyH2wLrBUppLxvjqIolcEvYWXsVTYpI6Hl1RWXH14Pq16xplsAx3dBo3YV12+cbH0e21l1MipKzSPMGerFcUwqU15lYWVKdXrh8lYBg572XO7a+X0KYhpBp9utnd4VEd3Qeua4KhtF8UCBiFQQbgHwPazMYLDCyhhjnGFljlbUqM3TqALPQa8iy+SP6vswdLyl4vrGwGWdoaig6mRUlEsPv4SVsXdk879NoPcDMPAp9+UlM2Oey7LOoxt4bnPsgqqTT1EuQUSkwBjjDCsTCrztDCsDrBWRhVhhZf5hjHkQa4rlNayMvcqmqKDijZhgKRsTYvnjrHkTvp4GvzsMEbUrfkZRlIvGH2FlbJtGuXYbeFI2zbrA77Og9SDLDhN3uXdF89+t8Nb1cGhdlcmqKMrFY5uyCTGOEZcnZRMSUrw14dhOzytRTsJrw8HVltJRFCVosG9k4zzx5NQnAv96AH6cB8d3+6Zs4i6HsEhd/laUIMM2m02RYBmHL+tUcSVjYOuncPqQFdWvUXvvDYeEQoM2lk+OoihBg43KxsC1z3iv6FyRGv4KJPb1rfGGbdRmoyhBhn2rUQbLfyYixvNmTKeySful720n9IKCnOKQFIqi2I5tyiYqVOCPifCLDyBpWMUVI+NgzxJrk2XjJN8a7363dbih9/SvSoWd0Hg3ihIYbFM2Pi19Q3Eqly8egzvLOjF6wY2P0aFT2WRMH+q6dgbZUhTFv9i39O08qWgHt5MbX7OmWr4Yh50UFsBLneA/f7xQ8RRFqWKC26kP4OdDkHfWt2VvJ6GOAduxHRckm6IoVY+NfjZOpz4PfjYAq161PuuUDafhhYbt4JgufytKsGCbzSa30MCAJ6CuFyWS0NuKvhfvaUe8Gxq1hX1f06JuRCm7TPO4qAuQVlGUi8UnZWOMGQy8hLUD9E0Rme6mzm1Y0bsE2CQid3hqM7fQQL9HvXeeNAyevoAwnw3bQmEuK+67worupyiKrXhVNr4EPzbGtAEeB3qLyEljTGNv7YYYgdOZENO02MZSlTTvamXNNKGeq2naF0UJCMZLCAqMMT2Bp0Xkesf14wAiMq1EnT8Bu0TkTV87bhwbIUcfjoIHt3mfSgWIxCmLSi2LK0pNxBhzXkSiA92vLwZiX4IftwXaGmO+M8asdky7POJajfK0EfNiKSosDrqlKIqt+DJ/8SX4cRjQBuiPFa90hTGmo4iUyp3rCKw8CaBJjGN642mrwsXy/s2Qdx4mLvGpeskplfNap1WKUjX4omx8CX6cCawWkXwg3RizE0v5rClZSURmAjMBLqsbYSksb0vfF0P9VrDlY8uTuGxucDeUVSzqXawoVYcv0yhfgh8vAAYAGGMaYk2r9nns2JfgWRfLZZ2thHUabkJRbMershGRAsAZ/Hg7VsrNrcaYZ4wxwx3V/g1kGWO2AcuB34qIR2PJz3khMHh6cSQ+f5DYx/pM/8Z/fSiK4hNeV6P8RUh4pBTl5/i3ExF4MRniu8Ft71b6cV2dUmoidq1G2eZBHB4qVqjPhm3814kxcN1zENPEf30oiuITtu2Nuqx2EcwKwKih482QWKmME4qi+AF7A577cyXKiYhlszm4xntdRVEAa4uSMWanMWaPMWZKBXVuM8ZsM8ZsNcZ86K1N26ZRIQbvsWyqAmNgwX3QrBP84n3/96co1Ry/bVHyl8DeMAb/LnuXpGVfyPhWYxIrim90B/aIyD4RyQPmAiPK1LkbeFVETgKIyFFvjdo8jQqQsknsawVXP6qJ6xTFB/yyRck2ZXM0OwT6PRaYzlo6UsCkrwhMf4oS3IQZY9aWOCaVKa/sFqVRwJvGmDiPnV6otBfLmTwD7YcEprO68VCvpZWWt+evAtOnogQvBSLiKRpdlW1RKoltI5uoMIHjewLX4dgFcMtbgetPUaov/tmi5AdBfSIxthCW/D5wHdZL9O8Oc0WpIfhri5K9eaMC+eMvLIAvn4TmaZB6W+D6VZRqiIh8Dnxe5t5TJc4FeMhx+IR9ygYC49TnJDQM9i6D7f+ybDgJvUqX714C2xdCQh9oNxgi6wZONkW5BLDXzyYQTn0lufHvltJ5Z4g1ynFuBF0yFT4YCZvmwqeT4M+t4V8PBFY2Ranh2OdBDIEd2QDEXwn3fmcpmpUz4GSG5VXcYQTENIZuE+HIJti6AGrXD6xsilLDsS3ERL2YCDm56wdo1tmW/tm91FIozdMqrNJ7+ldEnN5Hulj5xjVMqFITuORCTJzODbFP0QC0GeS1yneD9sOi38LkVdC4vYYJVZSLwDabTZ2IIjiRblf3vtFhBETEwLJn7JZEUao9timb1nULYfN8u7r3jdr1ofcDsHMRHFhttzSKUq2xTdkAgV+NuhB6TLYi/S2ZSvntIYqi+Iq9yiZQu74vhoho6D8Fju8k3hy3WxpFqbbYZiAGqoeyAejyS+h4C/K3dZrETlEuEFU2vhAaDqF1+e6xAVCY50oZrKtTiuI7tk2jdp8OhSuq0agg7xzMSIOVL9stiaJUS2xTNmfyQqBu2eBfQUxENMQ2g41zrCDqiqJUCtuUTVytIjjzX7u6vzA6j4ITeyFzrd2SKEq1wzZl06pOIfx3s13dXxhJwyEsCjZ5zVqhKEoZ1EBcGSLrQNL/wJaP4fpppYp6T/+KQ6eyXde6UqUopfFJ2Tgip78EhAJvisj0CuqNBOYB3UTE+1wj0Lu+q4Ke90GH4RASWur2oVPZpfKC60qVopTGq7LxJWGVo14s8P+A733vvZqNbMDaPGrnBlJFqab4MrJxJawCMMY4E1ZtK1PvWeBPwCM+917dplFOzvwX1r1DQy6vsErzuCh1AFSUEviibNwlrLqqZAVjTBeghYh8ZozxSdnsPBVmpVepjuScgq+ncVPoHcAdbquUVSw6rVIudXxRNh4TVhljQoAXgbu8NmQlw7ISYoWEQURtn4QMOhq1gxZXcfv+5ZbPjXH3ihRFKYkvS9/eElbFAh2Br40xGUAPYKExplwSLBGZKSJdRaRrw9oG8rPLVqk+pI3lipAjGnpCUXzEF2XjMWGViJwWkYYikigiicBqYLi31aiE2MLqrWySb+IcUcx/83kSpyyieVyU3RIpSlDjVdn4mLDqwgirhkvfTiKiie7+S0Z2bkrGtCFq/FVqFMaYwcaYncaYPcaYKR7qjTTGiLuZTFl88rPxlrCqzP3+vrQJVN/VKCc3/EntNUqNw1/uLvYGzwqx14H5onEqmlMH7JVDUaoWl7uLiOQBTneXsjjdXXJ8adQ2ZVMk1IxRwbaF8LcUOLzRbkkUpapw5+5SKkRDSXcXXxu1TdnsOFXNRzVOWl4NYZGwfrbdkiiKr4QZY9aWOCaVKffV3eXhSnVaeTmrhpyCGjCqAYiKgw43wuZ5MPD3EFXPbokUxRsFIuLJoFsZdxeApljuLh5Xoe212dQUet4HuWfg2xftlkRRqgK/uLuosqkKLkuFTrfDxg+rt++QouA/dxfbcn2HhEdKUb5PRuzqwdmjYEIguqHb4sQpi0qFoFAUu7jkcn3XOGIaW59FRZB7Wm03ilIGVTZVzYe3Wpszf/lJqdslQ05ouAnlUkSVTVXTagB8+QTsXQ5XDHDdLqlcNNyEcimiBuKqpvvdEHc5fPYby46jKAqgyqbqCasFt7xtKZr3b4Gcn+2WSFGCAlU2/qBFN7httpWqN+e03dIoSlCgNht/0eZay34TGmatUCHlMjIoyqWEjmz8iVPRLLgX3rsJTqTbLZGi2IaObPxNSAgk9IJ/Pwmv9YRrniCUBI+PeEp4p8nwlOqKKptAcOVd0PpaWPQwfPkki6KuYNDjR9gj8W6rN4+LqjDhnSbDU6orqmwCRd3mMGoObP2U9l89y9JxIyC2KZw9BrUbWCOgCijrEKgo1RFVNoHEGOh4MyTfVBw4bN6dcD4L+jwEHW+x7Dxl0GmSUhNQA7EdOBWNCHQdb23g/HQSzEiDtW9DQa698imKH1BlYyfGQMpIuPc7uP1Dazr12YOw9h27JVOUKkenUcFASAi0HwrthsC+5XB5T+v+7qXWZ+uBNSNes3JJoyObYMIYuOIaCHcYgVfNgA9ugQ9/oT46SrVHRzbBzB3z4IeZ8PU0eK0H9H0Yej/gc3I/9clRgglVNsFMWAT0ut9awfr372D5/0LDtpB8Y4WPlFQwnvx1FCXQqLKpDtRpBrfOgh73QbwjKP6uf0PT1FI+OFBewShKsOCTsjHGDAZeAkKBN0Vkepnyh4CJQAFwDBgvIvurWFalRTfrMz8bFjSOyN4AABN1SURBVPwKCnL5bsDj0P0et/45ihJMeDUQl8j7ewPQARhljOlQptoGoKuIpALzsVJyKv4iPAomLoHLe1jTq5n94IBP6ZYVxTZ8+e/QlfcXwBjjzPvrSjIuIstL1F8NjKlKIRU31G8Fo+fB9n/B4inwzmD49Trrvo+oAVkJJL4oG3d5f6/yUH8C8MXFCKX4iDHQYbi1XL77y2JFc3ANNL/S434r0E2dSmDxxc/GY97fUhWNGQN0Bf5cQfkkZ35hKSr0XUrFM7VirBUrgKPb4e3r4J0b4KctlWrGaWxOnLKI3tO/8oOgSnXBGDPYGLPTGLPHGDPFTflDxphtxpgfjTHLjDGe46bgm7LxlvfX2fkg4AmsNJxuN/eIyEwR6SoiXY1GrfMPDdvB8BlwfBe8cTV89pC1s9wHvptyDRnTh5IxfWip6ZVyaeEvO60vysZj3l+HcF2AN7AUjaYUsJOQEOgyxrLfdB0P62bB672gJmUfVfyNy04rInmA007rQkSWi8h5x+VqrEGIR7zabESkwBjjzPsbCrztzPsLrBWRhVjTphhgnrH28BwQkQvOCaxUAbXrw9C/wFX3wJFNEB4JCGz5xNqH5cUL2Z3/jhqPLxn8Yqf1yTlDRD4HPi9z76kS54N8aUexgYZtrAO40uyC+X+AOvFw9cOE06DCx8oqFjUe1yjCjDFrS1zPFJGZJa4vxE7bz2unlRJRqdask7Yw5hNrr9VnD7KiVj1YsQeumgwRte0WTwkcBSLS1UN5Ze20/Sqy05ZElc0lhbHCVVxxDez9il3vTqXp9zOh56+t4nNZEF3xaOdCUX+eaofLTgscwrLT3lGyQgk77WBf7bSqbC5FjKV0xubnkHFfL2vDZ1EhvNHXCuDVZQyk3GrZfaoA9eepXvjLTqvK5lInqp71WVQAvX8DG9+HLx6Ffz8BrfpZsZETe6vB+BLDH3ZaI+LW7uN3QsIjpajMcmx+fj6ZmZnk5OgyrT/IPJlNfL2oCq+d1Dq5m7r7FxObuZyf0h7hXLNeRJzeR+yhFZxrehU59dry0895FBS5/+6EhRia1o2sdL+BIDIykvj4eMLDw23pPxgwxpwXkeiA9xtMyiY9PZ3Y2FgaNGiA0TCYVc6OIz+TV1jkuo4IDaH9ZXUqfkDEOkJCYM2bVt4rgOhGkNgHWvSAtF9CROnv7Y+Zp0iNj/PputIyXQQiQlZWFmfOnKFly5Z+6aM6YJeyCappVE5ODomJiapo/ESlf8TGFMc+7jYR2v8P7P3KOvavhJ1fQNdxVvm6d+HMT3D5VYSYNkBchc2WJK+wqJwi8hfGGBo0aMCxY755VCtVS1ApG0AVTTAT2wQ6j7IOsLZBOJ0D96+EHz8ChGQTAo3aQ+tBcN2zVnlRkdeNoReLL6Mk/X7ZR9ApG6UaEdOo+PzmN+CGP0LmWrK2ryDi6EYKjh0hM/MUEaEh8HpPCK8NzdOoF9UOwvtAo3Zeu6jMNCuQoySl8qiyUaqOqDhoM4iGbYoXKuqDNappez0cWg+bPqJF3hn4BmvvVuepVvm2T+GyzkSENCylJCJCQypUIO4UkRK86L9OAOnVq9cFlXkiJiam0s+cOnWK1157zef6mzdvJiEhgddff73U/ezsbPr160dhYSGLFy+mXbt2tG7dmunTi6PG5uXlcXX//hQMeAru+gymHGDfbcs50O9F9jQdYimIrD0wfzzMSKP97FRSv7qL1J0vkxpxhPaX1WH8+PE0btyYjh07lurfOZJJjY+jnpxh8h0jSEpKIjk5mZdeeqnS70XxL0G1GrV9+3aSkpJc12U9Ty+WYPQNERFEhJALtGfExMRw9uzZSj2TkZHBsGHD2LLF93g3q1at4qGHHmLVqlWue6+++ioFBQXcf//9tG3bliVLlhAfH0+3bt2YM2cOHTpYUQn+8Ic/0Lp1a0aPHu2+8cICOLoNDq+3Rj+H18N/t8FtsyFpGBs+e5OWmR/z3rItDH1oBq06XQ2165da1Tpy5AhHjhwhLS2NM2fOcOWVV/LHv8/mpmt6lOtu+/btTPznkUvWq9mu1SjXlz3QhwmrJWXZtm1bqeuExz4rV+di8LW99957T7p16yadOnWSSZMmSUFBgaSnp0u7du1kwoQJkpycLHfccYcsWbJEevXqJa1bt5bvv//eVWfs2LGSkpIit9xyi5w7d87VbnR0tIiIpKenS/v27WXy5MnSuXNnycjIcJW9++67kpKSIqmpqTJmzBjXsyNGjJC0tDTp0KGDvPHGG+XaLMnZs2dlyJAhkpqaKsnJyTJ37txS5b/4xS8kMjJSOnXqJI888ohP72TPnj0SGxtb6l7Pnj0lPT1dVq5cKdddd53r/vPPPy/PP/+863rjxo1yww03+NSPi7zzIvk51vmuLyX3r51EptYpPl5Mke1b1jv+4OMiOT+Xenz48OHy9w8+cdv0tm3byn0Xek1bJgmPfeY6ek1bVjl5qxHAObHhN6/Kxo0Mw4YNk7y8PBERmTx5srz77ruSnp4uoaGh8uOPP0phYaGkpaXJuHHjpKioSBYsWCAjRoyQ9PR0AeTbb78VEZFx48bJn//8Z1fbJZWNMUZWrVpVqmzLli3Stm1bOXbsmIiIZGVlucqd5+fPn5fk5GQ5fvx4qTZLMn/+fJk4caLr+tSpU6XK09PTJTk52eu7KMnIkSMlIiJCMjIyREQkNzdXmjRpIiIi8+bNkwkTJrjqzp49W+677z7XdUFBgTRs2LBcm3369JFOnTqVO5YsWVKubnp6unTvlCR7vl8ksuJFkY/Gyo/pR6zCL38vMrWuyMtXiswbJ1kLfy+jezSVNTv2y6aDJ2XTwZOy/fBpV1vulE1Zqvq7F0zYpWzUQFyGZcuWsW7dOrp1s9KmZGdn07hxY66++mpatmxJSkoKAMnJyQwcOBBjDCkpKWRkZADQokULevfuDcCYMWN4+eWXeeSRR8r1k5CQQI8epYf4X331FSNHjqRhw4YA1K9fvDfp5Zdf5tNPPwXg4MGD7N69mwYN3G+aTElJ4ZFHHuGxxx5j2LBh9O3b9yLeCCxevJhz584xdOhQtm7dSkJCAsePHycuzprCiJupeMkl5tDQUCIiIjhz5gyxsbGu+ytWrKiUHOcKQjjXrBfED7H6dRqLk4ZbK10/babowA/U//lj3hoSQ622LSw/oe9eJuvYT+yvl0ROgw5kni6ieZw9HsyXMqpsyiAi3HnnnUybNq3U/YyMDGrVKg44FRIS4roOCQmhoKAAKO/HUZFfR3R0+SmziLit//XXX7N06VJWrVpF7dq16d+/v8ctHW3btmXdunV8/vnnPP7441x33XU89dRTFdb3RE5ODo8++igLFy7knXfeYcuWLQwZMoSoqCiXDPHx8Rw8WBxrKTMzk2bNmpVqJzc3l8jIyFL3+vbty5kzZ8r1+Ze//IVBg9xvvYkIDXGtSLlWn+K7QnxX8vPzGTZsGMOvG8V9o4cXOyQeWkuDHYtoUGT9G7UJjeTa5BsBh41m/0qIbgz1EjX/lh/RN1uGgQMHMmLECB588EEaN27MiRMn3P4gKuLAgQOsWrWKnj17MmfOHPr06VOpvm+66SYefPBBGjRowIkTJ6hfvz6nT5+mXr161K5dmx07drB69WqP7Rw+fJj69eszZswYYmJimDVrVqny2NhYn/+m5557jrFjx5KYmEhKSgoLF1oRYevVq0dhYSE5OTl069aN3bt3k56eTvPmzZk7dy4ffvihq42srCwaNWpUbj9SZUc2ULEXtIgwYcIEkpKSuO/hJ0oX3jYbCnLh2E74aTOntv2H+o6AYgDMHQ3ZJyA0AupfAY3aMi62JYmOMN8t64ay/PHBlZa1stT4UBx2zN0kiG02IiJz586VTp06SUpKiqSlpcmqVavK2TnuvPNOmTdvnogU20DS09MlKSlJ7rnnHklJSZGbb765QgNxWZuJs2zWrFmSnJwsqampcuedd4qISE5OjgwePFhSUlJk5MiR0q9fP1m+fHmp50qyePFiSUlJkU6dOknXrl1lzZo1IiJyww03yKFDh0REZNSoUZKcnOwyEJcsc7Jjxw7p3r275Ofnu667dOniKh8/frzLvrJo0SJp06aNtGrVSp577rlS7cybN08eeughj+/cE7fffrs0bdpUwsLCpHnz5vLmm2+Wk3nFihUCuP7uTp06yaJFi9y2V+p7VlQkcnCNyPr3LdvPB78QeamzyJdPWeW556RgapwceSpBvv99N5n/5FD5x7OTRDLXFj9fVHTBf1tJyn4//WU3wiabjS59VyEXsqRcndmwYQN//etfee+99zzWu/nmm5k2bRrt2nn3GA4EZb9nbhGxpmE5p+H7N+DkfjiZASfTKTp9mJAbpkGPyXB8Dzmv9GR/UWMOSGMypRGHpCHLitJIl8sIxUpZVIiVTcTTdzBxyqJycX/8kbddN2K6oUYNIWsgXbp0YcCAARQWFhIa6j41T15eHjfeeGPQKBqfcdp7IutCv0dLFQ2ctpifFpwje8EimnGc/xd9Pbe3K6LdyQw4+S3kn+PJ266HDkNh39fw3s0Q2xTqNGfZ4XBmPvkmHxQOZL80pQ7nqGPOcUziaB5XN+B/ZiAJ6pGNovgDv37PRKzRUFgtKyf78d2waS78fAhOZ8LPh63zuz6H+Cthwwfwz19Zz9aqAzGNIaYJjHiVxD9tI+OBFvDTj9Y9Z1nthuUM2d7sPSXL9/9xmI5sFKXaY4y1R8xJwzYw8Pel64gjThBAQk8Y8Sqc/S+cPVr8GV6b5nFRvPDKSzwcPr9sJ/DIbmsj7Pr3YNdifn32LLcP7ARR9aF2A1r9n0OGc8cB+OnUWTKmW1E7zR/98Hf7gCobRQk0JeME1W9VnKO9DN9NaQL5PeHs7+HsUR6fvZTQ80dpZE7x0nOrKSKECaGruTV0E4PCzsLqlVCYByFhFPGu1cjSqbDhffZGAtPrWjGmbUKnUcolR439nolA3jnIPknvv+/m0KlsupkdJIXsJyEyhwlX1oXzJzC3vq3TKKjYsU1RqgK7/nMNCMZArRioFcN3U5xpn9ytZr0dSKlcBFWIicjISLKysmr2F0KxDRErBnFZT2YlMATVNEqzKyj+RrMrBHl2BWPMYOAlrIRVb4rI9DLltYDZwJVAFvALEcnw1KY7ZaMoiv/xRdn45Tfvg2ChwKvADUAHYJQxpkOZahOAkyLSGngRsGlxTVGUi8Vfv3lfbDbdgT0isk9E8oC5wIgydUaAc62N+cBAo1ZeRamu+OU374uyaQ4cLHGd6bjnto6IFACnAfsW9BVFuRj88pv3ZenbnbYqa+jxpQ7GmEnApBLX533oP5CEAQV2C+GGYJRLZfKNYJSptjFmbYnrmSIys8R1lf3mS+KLsskEWpS4jgcOV1An0xgTBtQFTpSTxPqDZgIYY9aKSFcf+g8YwSgTBKdcKpNvVFOZquw3XxJfplFrgDbGmJbGmAjgdmBhmToLgTsd5yOBr0SdZRSluuKX37zXkY2IFBhj7gf+jbUM9raIbDXGPAOsFZGFwFvAe8aYPVja7fZK/GGKogQR/vrN+7RdQUQ+Bz4vc++pEuc5wK2+/jEOZnqvEnCCUSYITrlUJt+oljL54zdvmwexoiiXFkG1N0pRlJqL35WNMWawMWanMWaPMWaKm/JaxpiPHOXfG2MSg0Cmu4wxx4wxGx3HxADI9LYx5qgxxm0AY2PxskPmH40xaUEgU39jzOkS7+nC8sVUTqYWxpjlxpjtxpitxpgH3NQJ6LvyUaaAvitjTKQx5gdjzCaHTH9wUyewvz1/RlPHMi7tBVoBEcAmoEOZOr8C/u44vx34KAhkugt4xZ9yuJHraiAN2FJB+RDgCyz/hh7A90EgU3/gswC/p8uANMd5LLDLzb9fQN+VjzIF9F05/vYYx3k48D3Qo0ydgP72/D2yCcatDr7IFHBE5Bs8+ymMAGaLxWogzhhzmc0yBRwROSIi6x3nZ4DtlPduDei78lGmgOL42886LsMdR1kDbUB/e/5WNsG41cEXmQBucQzB5xtjWrgpDzS+yh1oejqG6l8YY5ID2bFj2N8F63/tktj2rjzIBAF+V8aYUGPMRuAosEREKnxPgfjt+VvZ+MXt+SLxpb9/AYkikgospVj720mg35MvrAcSRKQTMANYEKiOjTExwMfAb0Tk57LFbh7x+7vyIlPA35WIFIpIZywP4O7GmI5lRXb3mL/k8beyqYzbM766PftbJhHJEpFcx+U/sGJ22I0v7zKgiMjPzqG6WH4Z4caYhv7u1xgTjvWj/kBEPnFTJeDvyptMdr0rR3+ngK+BsjmEA/rb87eyCcatDl5lKjO/H441B7ebhcBYx0pLD+C0iByxUyBjTFPnHN8Y0x3r+5Tl5z4NlvfqdhH5awXVAvqufJEp0O/KGNPIGBPnOI8CBgE7ylQL7G8vAFbxIVjW+b3AE457zwDDHeeRwDxgD/AD0CoIZJoGbMVaqVoOtA+ATHOAI0A+1v84E4B7gXuleHXhVYfMm4GuQSDT/SXe02qgVwBk6oM11P8R2Og4htj5rnyUKaDvCkgFNjhk2gI85eZ7HtDfnnoQK4oSENSDWFGUgKDKRlGUgKDKRlGUgKDKRlGUgKDKRlGUgKDKRlGUgKDKRlGUgKDKRlGUgPD/AYOoW2nVAVHkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.hist(\n",
    "    x,\n",
    "    bins,\n",
    "    density=True,\n",
    "    fill=False,\n",
    "    histtype=\"step\",\n",
    "    label=f\"empirical s.t. $\\lambda(0)={lamb(0)}$\",\n",
    ")\n",
    "pl.legend(loc=\"lower left\")\n",
    "pl.xlim(xmin=0, xmax=xmax)\n",
    "pl.twinx()\n",
    "\n",
    "pl.plot(\n",
    "    bins,\n",
    "    mixture_pdf(bins[:, None]),\n",
    "    \"--\",\n",
    "    color=\"C1\",\n",
    "    label=f\"mixture pdf\\n$\\hat\\lambda(0)={lamb0:.3f}$\",\n",
    ")\n",
    "pl.legend(loc=\"upper right\")\n",
    "pl.ylim(ymin=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0.28398601412773133)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASEAAAD4CAYAAACjW1BIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deXgUVda435uEEAjIJiAQIWwCZgWDgKwOiArKoqAgERQQ12+ccXDEz++nDuMoozgy7iJuyCrICCgwyqKCgOyC7EIiBDcIq4QsnZzfH9VpOkkvlaTTVUnf93nqSVXd7XSl6/Rdzj1HiQgajUZjFWFWC6DRaEIbrYQ0Go2laCWk0WgsRSshjUZjKVoJaTQaS4mwquGwsDCpUaOGVc1rNCFLVlaWiIhtOiCWKaEaNWpw/vx5q5rXaEIWpdQFq2VwxzbaUKPRhCZaCWk0GkvRSkij0ViKZXNCmqpNXl4eGRkZZGdnWy1KyBIVFUVMTAzVqlWzWhSfaCWkqRAyMjKoXbs2sbGxKKWsFifkEBEyMzPJyMigZcuWVovjEz0c01QI2dnZNGjQQCsgi1BK0aBBg0rRE9VKSFNhaAVkLZXl+WslpNFoLEUrocpAfp7VElRZlixZwpQpU3zmef/99/npp5+CJJFnvvzyS2666SYAcnJy6NevH8nJycyfP99SuQKBnpi2MwX58PYfoGVP6P8MAN2nrObY6YsGr83q1uCbSX+wSsJKz6BBgxg0aJDPPO+//z7x8fE0bdrUdL0Oh4OIiIp5vbZv305eXh47duyokPqDje4J2Zlzv8DPO2DLe1BQAMCx0xdInzLQdbgrJM1F0tPTad++PePHjyc+Pp5Ro0axcuVKunfvTtu2bdm0aRNgKJiHHnoIgMGDBzNz5kwA3nrrLUaNGsXChQvZsmULo0aNIjk5mQsXLhAbG8uJEycA2LJlC3369AHg6aefZsKECfTv35/Ro0eTn5/Po48+SufOnUlMTOStt97yKueYMWNITExk2LBhZGVlAbBixQrat29Pjx49WLRoEQC//fYbqamp7Nixg+TkZA4dOlShzzEYaCVkZzI2G39zf4cj662VpRLyww8/8PDDD7Nz50727dvHnDlzWLduHVOnTuXZZ58tkX/69OlMnjyZtWvX8uKLL/LKK68wbNgwUlJSmD17Njt27MDfpuutW7eyePFi5syZwzvvvEOdOnXYvHkzmzdv5u233yYtLa1Emf379zNhwgR27tzJJZdcwuuvv052djb33HMPS5cuZe3atfzyyy8ANGrUiBkzZtCzZ0927NhB69atA/OwLEQrITtT4Lh4vmuhdXJUUlq2bElCQgJhYWHExcXRt29flFIkJCSQnp5eIn/jxo2ZPHky1157LS+++CL169cvdZuDBg1yKarPP/+cmTNnkpycTJcuXcjMzOTgwYMlylx++eV0794dgNTUVNatW8e+ffto2bIlbdu2RSlFampqqWWpLOg5ITtTqISaXQV7PoEbn7dWnkpG9erVXedhYWGu67CwMBwOh8cyu3btokGDBj4noiMiIihwDo+L2+FER0e7zkWEV155heuvv96nnMWX0guvK8sSe3nRPSE7U7gq1u0huGEKoCOjVCSbNm1i+fLlbN++nalTp7qGTrVr1+bcuXOufLGxsWzduhWAjz/+2Gt9119/PW+88QZ5ecb/8cCBAx7d1xw5coQNGzYAMHfuXHr06EH79u1JS0tzzfnMnTs3MB/ShmglZGcKnEqoeTdIGgER1X3n15SZnJwc7rnnHt59912aNm3Kiy++yNixYxER7rrrLu677z7XxPRTTz3Fww8/TM+ePQkPD/da5/jx47nyyivp1KkT8fHx3HvvvR57YB06dOCDDz4gMTGRkydPcv/99xMVFcX06dMZOHAgPXr0oEWLFhX58S1FWRV3LDo6WrRTMz/8fhxO/wiXJUL2Gdg5nw5LmrB3yq2uLLGTPiN9ykALhfTM3r176dChg9Vi2J709HRuuukmvv/++wqp39P/QSmVJSLRXooEHd0TsjO1GkJMCkREwon98PkTXBe2zWqpNJqAoiem7cwvu3huxlze+b0L+YSxvnp9bovaaLVUmgASGxtbYb2gyoJWQnbm4Oc87niNx//+JFSLgv9upMm3b8KFU1CjntXSaTQBwdRwTCl1g1Jqv1LqB6XUJA/pLymldjiPA0qp04EXNQQpyDf+hjudUsXfYizb71tmnUwaTYDxq4SUUuHAa8CNwJXASKXUle55ROTPIpIsIsnAK8CiihA25ChcolfOf1PTTlCvJWSWNHjTaAKJiY7HI0qpPUqpnUqpVUqpFm5p+W6dkiX+2jIzHLsa+EFEDjsbmAcMBvZ4yT8SeMpEvRp/FDjIlXAiC43WlIIHNhpDs8rGex5W8OKGwNX3QG4WzB5eMj35Dug4Cs5nwkeji6bd/ZmpZnft2sV1113HypUriY+PL4PgoYdbx+M6IAPYrJRaIiLu7/x2IEVEspRS9wPPA7c70y44OySmMDMcawYcdbvOcN7zJHwLoCWw2kv6BKXUFqXUFm8Wqxo3CvLIp5gdSqECKhyqaXzy7LPPsn79eo97xTRecXU8RCQXKOx4uBCRNSKS5bzcCMSUtTEzPSFPtuPejItGAAtFxOMbIiLTgelg2AmZkjCUueZhblrTnFXF7y9+EM6fgDsqkS8ZXz2XyJq+06MbmO75FKfQ0njOnDllKl9FiVBKbXG7nu58Nwvx1PHo4qO+ccByt+soZ/0OYIqIfOJTGBMCZwCXu13HAN421owAHjRRp8YMtRpySDx0Oms2gO/mQdbJ4MukqQo4RCTFR7rpjodSKhVIAXq73W4uIj8ppVoBq5VSu0TEq88RM8OxzUBbpVRLpVQkhqIpMdmklGoH1AM2mKhTY4ZDqxkR7mFkGzfUuUpWtt5BKLFr1y5atGjBG2+8UeT+hQsX6N27N/n5Rqd97NixNGrUqMi8UW5uLr169fK62dUMK1asoF27drRp08anB8d///vfxMfHExcXx7Rp00yVf+mll4iLiyM+Pp6RI0cG0qm9qY6HUqof8AQwSERyCu+LyE/Ov4eBL4GOPlsTEb8HMAA4ABwCnnDem+xsvDDP0xhdL1N11qxZUzR++OQBOfZkbMn7BQUiLyWIfHiLtHjs0+DLZYI9e/ZYLYKL9evXS9euXYvce/XVV2XatGmu66+++kq2bt0qcXFxRfI9/fTTMmvWrDK163A4pFWrVnLo0CHJycmRxMRE2b17d4l8u3btkri4ODl//rzk5eVJ37595cCBAz7LZ2RkSGxsrGRlZYmIyPDhw+W9994rUben/wNwXny/7xHAYYz53UjgOyCuWJ6OTn3Qttj9ekB15/mlwEHgSl/tmbITEpFlInKFiLQWkX847z0pIkvc8jwtIiWW8jTlIN+BQzxskFTK6A0d/pK6nCuZrilCo0aN2L17d5F7s2fPZvDgi3OtvXr18ug/aMiQIcyePbtM7W7atIk2bdrQqlUrIiMjGTFiBIsXLy6Rb+/evXTt2pWaNWsSERFB7969+c9//uO3vMPh4MKFCzgcDrKyskrlftYXIuIAHgL+C+wFPhKR3UqpyUqpQl+4LwC1gAXFluI7AFuUUt8BazA6Jt5W0gFtMW1vCvJwFF8dKyRpBNSsT8HS0PA5Ux4mTZpETk4OP/74Iy1atCA3N5fDhw8TGxvrt2x8fDybN28ucb9nz55F3HsUMnXqVPr16wfAsWPHuPzyi6OamJgYvv32W49tPPHEE2RmZlKjRg2WLVtGSkqKz/LNmjVj4sSJNG/enBo1atC/f3/69+/v9/OYRUSWAcuK3XvS7byfl3LrgYTStKWVkJ0pcHhXQo06QKMOnF2q54V8sWLFCs6fP8/AgQPZvXs3LVq04MSJE9StW9dU+fDwcCIjIzl37hy1a9d23V+7dq3fsuLBQ4UnR2UdOnTgscce47rrrqNWrVokJSURERHhs/ypU6dYvHgxaWlp1K1bl+HDhzNr1qxK6YFR76K3M/mOknZC7uT8zqCwb/QqmReys7P561//yuuvv05CQoJro2iNGjVKNYmbk5NDVFRRA9GePXuSnJxc4li5cqUrT0xMDEePXlzpzsjI8DpkGjduHNu2bePrr7+mfv36tG3b1mf5lStX0rJlSxo2bEi1atW45ZZbWL++cvoh1z0hOzP4VUZ/t4ySgwEnmT/wcuRrsLcdXDUmmJJVCp555hlGjx5NbGwsCQkJLFliTFvUq1eP/Px8srOzSyiX4mRmZrpedHfM9IQ6d+7MwYMHSUtLo1mzZsybN8+rvdJvv/1Go0aNOHLkCIsWLWLDhg3Url3ba/nmzZuzceNGsrKyqFGjBqtWrSIlxdequ33RPSE7U7M+x/GxW75JEj8WNDL8T2uKsH//fr744gv+9Kc/ARTpCQH079+fdevWua5HjhxJt27d2L9/PzExMbzzzjsArFmzhgEDBpRJhoiICF599VWuv/56OnTowG233UZcXJwrfcCAAS5f1rfeeitXXnklN998M6+99hr16tXzWb5Lly4MGzaMTp06kZCQQEFBARMmTCiTnJbja+msIg+9RG+C7bPlz4//1WeW155IFXm6nsjvJ4IklDnstETviW3btklqaqrffEOHDpV9+/YFQaKKoSxL9ME+dE/IzmybybDwr31m+Sy/K0g+7Ps0SEJVDTp27Mi1117rMlb0RG5uLkOGDKFdu3ZBlCz00ErIzvhaHXOyW1oY7j0yvM4cabwwduxYn47qIyMjGT16tNd0TWDQE9N2Jt+HnZALBfes1p4WNZUW3ROyMyZ6QgDUrG9YUWs0lRCthOyMWSUEsGoyzL+zYuUpJWJROCmNQWV5/loJ2Zl7VvOXvPvM59/3meFnyAZERUWRmZlZaV6EqoaIkJmZ6dcOyg7oOSE7ExnNBUx+ieKGwtoXYe8SSBlbsXKZICYmhoyMDI4fP261KCFLVFQUMTFldngYNLQSsjNfPc+NYb8DJiKsNo6HBm1h18e2UELVqlWjZcuWVouhqQRoJWRnNr1NzzCTztmVgoRh8OUUOPsTXNKU7lNWc+z0BVeWZnVr8M2kP1SQsBpN2dBKyM6UZmIaIGE45Fx0L3Hs9IUicepjJ+kd9xr7oZWQnSmtEmrQmu7b+3JszXZgO83q1qgw0TSaQKGVkJ0prRICfj59nvQJl0CDNlDHY2QmjcZWBCQMtDPPbc6IjLuVUjq+SiAogxJqxCmYORh2lM0lqUYTbAISBlop1RZ4HOguInHAnypA1tDj/37jBcdtpSryCw2gxTWwayFoGx1NJcBMT8hvNEbgHuA1ETkFICK/BVbMEEUpyuToIP5WOLEffv3ef16NxmICFQb6CuAKpdQ3SqmNSqkbAiVgyJKfB0sfplfYd6Uve+UQCIswekMajc0JVBjoCKAt0AcjUNpapVS8iJwuUpFSE4AJYLhJ0PjAkQNb36e9GukzW7O6NYosvTerW8MIm9zqWkjz7YtIo7EDgQoDnQFsFJE8IE0ptR9DKRVxciM6Fr15CvIAcPj5F3k1PhzyOtQoGUdLo7EbgQoD/QlwLYBS6lKM4dnhQAoachQYHv/ySrk65qJWIwiP0JPTGtvjVwmJuWiM/wUylVJ7MKIuPioimRUldEiQb/SEfIb88cfepfB6V8g9HyChNJrAY8pYUfxHYxTgEeehCQSSD9Wiycmt5j+vN2rUh+P7YO+nkHR74GTTaAKI9idkV+rEwBM/8XFBr7LX0bwb1IvVhosaW6OVUFUmLAyS7jBWyU4f9Z9fo7EArYTsyumjsHAcSeqH8tWTNAIQ2DkvIGJpNIFGKyG7cuEkfL+QxupU+eqp1wL6/C80vyYwcmk0AUYrIbtS4AAgLxCODvo8BrHdy1+PJmTwt2ldKfWIc8P6TqXUKqVUC7e0MUqpg85jjL+2tBKyK/mGEsoP1L/oxEFuCNsUmLo0VRozm9aB7UCKiCQCC4HnnWXrA08BXTD2nT6llPIZFE8rIbsSyJ4QwLppTK32JmSfDUx9mqqM303rIrJGRLKclxsxdlIAXA98ISInnRvavwB87iXVSsiuhIVDrcbkSDnshNzpPJZaKht2zg9MfZrKTIRSaovbMaFYuplN6+6MA5aXsaz2rGhbmneFiQfYFii/0M2u4ruCViRtngGdx4NS2hF+6OIQkRQf6WY2rRsZlUoFUoDepS1biFZCIcSs/H4kHZ8OP34DsT20I3yNN8xsWkcp1Q94AugtIjluZfsUK/ulr8b0cMyuHPkW5txOjAqcf7il+d2g1mXGVg6Nxjt+N60rpToCbwGDijkx/C/QXylVzzkh3d95zyu6J2RXzh6DAyuoQeCGR9lUhz/tggjty0njHRFxKKUKN62HA+8WbloHtojIEuAFoBawQCkFcEREBonISaXU37noxmeyiJz01Z5WQnbF6cojYEv0hRQqoPPayYHGOyY2rffzUfZd4F2zbenhmF1xOjUrsz8hX6z4X3izOxE4Al+3RlNKtBKyK05/Qg6pgM5qbA849zP9wrYFvm6NppRoJWRXIqOhfityK2LEfMX1ULc590R8pj0vaixHKyG7kjAM/ridTOoEvu6wcLjmj1wVdtBYrndS6DQ/dtJndJ+yOvDtajQe0EooVOmYynG5BDZNd936ZtIfSJ8ykPQpA4sYMWo0FYlWQnZl5wJ4/yaqk1sx9VerwbjcR2HIGxVTv0ZjkoDEoldK3aWUOq6U2uE8xgde1BDj9I+QvpaCCvyd2CmtjbknPS+ksZCAxKJ3Ml9Ekp3HjADLGXo4d9E7KrqzemybEZEj81DFtqPReCFQseg1gabAASqsbLHoS0OdGDiVDt9Mq9h2NBovBCoWPcCtTi9rC5VSl3tIRyk1odB9gMOhDeV8kp9nxJOvaGo1go6psGMunDlW8e1pNMUwo4TMbM1fCsQ6vaytBD7wVJGITBeRFBFJiYjQO0Z8UvsyaJIcnLa6PwxKwZfPBac9jcYNM0rI77Z+Ecl028r/NnBVYMQLYbreD+O/CE5bdZtD53uM+GTHDwSnTY3GiZnuiGtbP3AMY1v/He4ZlFJNRORn5+UgjHDRmspEr4nQtCM0aANcNFwsRDs801QUfpWQyW39f3TGpXcAJ4G7KlDm0OCr5+Hot8DY4LRXsz4kDjfORUooHO3wTFNRBCoW/ePA44EVLcQ5lQ6/WeB8bNtM2LUA7lxsRHDVaCoY/S2zK/l5EG7B5H14dSNs9PcfB79tTUiilZBdKXBAWIAibZSGhOFwWSKs+hvkng9++5qQQyshu1IQJDuh4oSFwY3/hDNHYc2zwW9fE3JoJWRXLm0HzSyydGhxDVx1N2x8A04f9Z9foykH2mLQrvT9f8bfjRatSl33N2NoVtej8btGEzB0T0jjmag6ENvdOP/9uLWyaKo0WgnZlY/vgUXFo/NawHfz4N+JtFIlYt9pNAFBKyG7cvoInPvZf76KptW1EB7J1GpvupzvazSBRCshu2LVEn1xajeGm6fRKewHWDXZamk0VRCthOyKVUv0nogbyoeOfrD+ZTjgM6KvRlNqtBKyKwX5EG6DnpCTZxyp0CzFHkNETZXCJj+1mhI07waXNIEdVgtikEMkjPvcCBek0QQQrYTsysCpxt/PbLR7vVAB7VkMP22Hfk+bKtZ9yuoiIYS0WxCNO1oJaUrP0U2w4VWo3QS63Os3+7HTF0ifMtB1rd2CaNzRc0J25c0e8N8nrJbCM9dNhnYDYcUkPVFdRTER5quXUmqbUsqhlBpWLC3fLfzXEn9taSVkV879at9d7GHhcOvbcFkCLBwLv+yyWiJNADEZ5usIhvPCOR6quOAW/muQv/a0ErIrdlqi90RkNIycD9Uvgb2fWi2NJrD4DfMlIukishMoKG9jNv6Whzj5Dlst0XvkkiZw31qo2cBqSTSlI0IptcXterqITHe79hTmq0sp6o9y1u8ApojIJ74yByQMtFu+YUopUUqllEJgjSfs3hMqJPpSI1zQb3th9m2QfcZqiTT+cRSG3nIe04ulmwnz5YvmIpKCERBjmlKqta/Mfr/lbuPD6zA04mal1BIR2VMsX23gj8C3pRBW4424odAkyWopXPiNvnH2GBxaDXNuh9RFEFnTAik1AcJvmC9fiMhPzr+HlVJfAh0Br3HGzfzUusaHAEqpwvHhnmL5/g48D0w0K6zGB0PfdJ7YYznbb/SNNv2MyeqFY2H+KBgxF6pFBVFCTQDxG+bLG0qpekCWiOQopS4FumPoBa8EJAy0UqojcLmI6BnKUCZuKAx6BQ6tgXl3QF621RJpyoCIOIDCMF97gY8Kw3w5Q3uhlOqslMoAhgNvKaV2O4t3ALYopb4D1mDMCRXvsBTBTE/I5/hQKRUGvISJWGNKqQnABIDIyEgTTYcoeRfg+Vbwh/8DYq2WpnR0TAURI2yQlHvhRGMRJsJ8bcYYphUvtx5IKE1bZpSQv/FhbSAe+FIpBXAZsEQpNUhE3GfgcU6ATQeIjo4uzURXaFHggLysyvsSd7oTkkcZTvOzz1KdXKsl0tiYcoeBFpEzwKWF186JqInFFZCmFBQ6D7ODPyEvFJ+oLp72zaQ/GJ4AZg9nRrXfIfd6PVmt8UigwkBrAklBvvHXxjvWfW1AdSmnsHC4agzdjzwAc0fAyHlaEWlKYMpOSESWicgVItJaRP7hvPekJwUkIn10L6icFDh7QnY3VjRD8h1MzLvXiOo693bIzbJaIo3NqATWcCFIRJQR96theyAzYNV6svUJBosKevGvoZ3gk/vg0z9RbAeAJsTRSsiO1KwPN09zXgTOTshSHz5Jt0NEpBFietM+123ta0ijlZAdETGOsCq2vzhuqPNkL2x8EzrdqX0NafQuelvy626YXA/2VM05/wSVBv99HGYPpybaoDHU0UrIjhROTFeGDaxlYJe0gltnwJGNvBf5POT8brVIGgvRSsiOuJboq6YSAiD+Vrh1BlepAzB7GOScs1oijUVoJWRHCo0Vw6uwEgKIv4U/5j0EP38HP++0WhqNRWglZEcKHMbfqtwTcrKsoCs8vBNiuwOgyu+oT1PJ0ErIjlzSFK75H6hzuf+8VYFaDY2/Oz/io8jJ2jFaiFH1f2orIw1aQ/9nnBc+vSBULarVIFkdgg9vgTsXQVSdUhXXNkeVE62E7IgjFxzZhjP5UKLDzTyY90em//wKfDgUUhfR/d/bXIrFn1LRNkeVEz0csyMHVsCUyw2/zSHG5wWd4baZxkT1h0M5cfoM6VMGkj5lYJFejqbqoHtCdqSS2wmVe49a+4Fw+4fw83fkHK4Cm3g1Pqmc3/KqTr5zdayS7qIPyDxMuxuNY8VnRo+w9mXlr1NjS7QSsiOVvCcUSKqTa0xU12pIHR60WhxNBaC/5XYkhOyEilNyKFcHBr0M80YxO/JZyLrO8DKgqTKE3re8MtAkCfo8DlGXWC1J0PE6lBsxh7azRsAHN8PoxUbQRU2VQK+O2ZGmHaHPJKhe22pJ7EPbfozLmwiZP8BX/7RaGk0A0T0hO5J9FnJ/h1p6MtaddQUJcPdyaNQB8GycqKl8mFJCSqkbgH9jOLqfISJTiqXfBzwI5AO/AxP8BTzT+GDzDFj1N3jiF6slsR/NOhl/L5xm0vl/cvPj70KdZr7LaGyN3+GYWyz6G4ErgZFKqSuLZZsjIgkikowR8vVfAZc0lAgFVx7l5VQafcK+g/cHwOmj/vNrbIuZOSFXLHoRyQUKY9G7EJGzbpfRuEVo1ZQBvUTvn6YdSc19HLJOGYro1I9WS6QpIwGJRQ+glHpQKXUIoyf0R08VKaUmKKW2KKW2OByOssgbGhQ4QIWD8hSBu/JTuAwfO+mzcs3jfCdtYMxiYw7tvQFw8nAApdQEi3LHonfdEHkNeE0pdQfwf8AYD3l0GGgz5OdVWmtpMwR0Z3vTjjBmKSx5CMIjA1evJmgEIhZ9ceYBb5RHqJCn3Y1wiZ5sNU2TRJjwldFzLMiHEwdcK2ga+1PuWPQASqm2InLQeTkQOIim7LS4xjg05ikcun49Fdb9y9iJ7wPte8g3JlbEewHTgERghIgsdEsbgzEaAnhGRD7w1VagYtE/pJTqB+QBp/AwFNOUgrM/Q16W4dxMUzo6j4P9y2DuSIaGTcD4TSyJ9j3kHbcV8eswRkKblVJLipndHAHuAiYWK1sfeApIwZi22eose8pbe6aWX0RkGbCs2L0n3c4fNlOPxiRr/gE/rIS/7POfV1OU6Evhrk9h3h28lPYGbGgO3R6wWqrKhmtFHEApVbgi7lJCIpLuTCvuFPx64AsROelM/wK4AZjrrTG9bcOOFDggrOpOTFc41WvDHQtYnt/ZMPo8fcRvEfcVu+5TVgdBSEuJKFyldh4TiqWbWhH3QqnLakMUO1LgqPrhfiqaalE8mPcwhx9oAXWbG/dEvJo9uM8HhcDQzCEiKT7STa2IB6qs7gnZkfw8bagYAAoIM1bOALbNhDm36SCL5ijtini5ymolZEf0cCwguA+xpizbAz+sgndvhDPHrBbN7rhWxJVSkRgr4ktMlv0v0F8pVU8pVQ/o77znFf1za0c6j9O/2AGg6BALJo3vBx/dBW//gTj1P9YJZnPMrIgrpToD/wHqATcrpf4mInEiclIp9XcMRQYwuXCS2htaCdmR1tpepUJo0w/G/Rfm3M6CyMlw9hYj0KSmBCZWxDdjDLU8lX0XeNdsW3o4Zkd+2wfHD1gtRdWkcRzcs5qnHaO1ArIJWgnZkc8egU//bLUUVZdajfgo/1rj/McNMG+UDj1tIXo4Zkfy8yCyptVSVCm8xkI7lWYEm5zex9jqcVmCNQKGMFoJ2ZECh16iDzBe94Ul3wH1WsLCu2FGPxgwFdDRPIKJHo7ZkYI8vUQfTFp0g3vXwuVdYMlD9A3barVEIYVWQnakIB/Cwq2WIrSo1RDu/A8MeZPVBR2Ne3kXfJfRBASthOzIdX+Hbg9ZLUXoERYOySMRwuDsT/ByR9j4hrHdQ1NhaCVkR9r2M4YIGusIrw5NkmHFJJg93HCvoqkQ9OynHflxPdS8FBpeYbUktsLrCldFEN0ARs41wi99/v/g9a7GpHXi8IprM0TRSsiOLBwLbWBRt6kAABTKSURBVK+DQa9YLYmtCLrnQ6Xg6nug1bXwyf3GUr5WQgFHKyE7onfR24tL28DYFRcnqn/dA7/shMTbq2xElGCiv+l2RC/RW4qnYd83k/4A1WsZNzbPgC3vwM75cNNLUC/WGkGrCFoJ2ZF8R5UO+WN3ig/7Sjg5G/ACNGxveG18vRv0mmisZkZUL1e7oep8P1Cx6B8BxgMO4DgwVkR0SMyyUuDQdkJ2JiwcukyA9gNg+WOwarIR8+ya8rkHCVXn+36VkEnP+9uBFBHJUkrdjxGF9faKEDgkGDkH6jS3WgqNE6/DszoxMGI2HFpjWFsDHPkWajYw5pE0pjDTEzLjeX+NW/6NQGoghQw52vSzWgKNG36HZ62dO/JFYNlE+G0vdL0fej0KUZcEScrKS8Bi0bsxDljuKUHHojdBQT7sXQqZh6yWRFNalIJRC41Vs/Uvw6spsH228T/VeMWMEjLtPV8plYoR9OwFT+kiMl1EUkQkJSJCz4l7xJED81MNRaSpfNRuDENeg/GrjOHa4geMYIwarwQsFr0zAusTQG8RyQmMeCFIQZ7xV9sJVW5iUmDcSiY+N5WPPwDhMwaFfYOKbsi/n9AO69wJVCz6jsBbwA0i8lvApQwlCrvueom+8hMWxsJzccaKlwi89Rz8sgs+WA19nzIUlSZgsehfAGoBC5RhQXpERAZVoNxVl/zCnpBeoq8sFLfvcce1v00pGL+Kvz31F576dTnM6AvtBkK/p6BhuyBKaz8CFYteL+cEigLnhL22mLYtnpbs3e17vBJRnffyb+Sph6fAt2/AN69w39srWXH2B6qTS7M6obmSpice7Eb0pXD3CqjfympJNF4otxVz9VrG8n2X+1nx1FeGAlv6J/h1N+zPhStuCKk9aVoJ2Y2I6tqXUKhQuBcNoGlHOLQK5o6ARnHQ48+EU75tIJUFrYTsxoVTsH8FxPaAupf7z6+pVHj1iXTVGMPp/veLYN2/YNF4/hIxCKj6U6taCdmNMxnwyX1w24daCVVBfA7lwqtB0u2QMBwOLOejD37mATC2gmyfCVffC00SgyVq0NDuXe1G4eqYXqIPXcLCoP1A0qWJcX1iv9FDeqsnvHsj7P7Pxe9JFUD3hOxGoZ2QXh3TFNJpNHS4GbbPgk1vw4K7oEkSTPjKNYHtzw2ILzMCq9FKyG4UaDshjQdq1DNchXR9AA5+DlmZhgIqyIeFY+l69jJenPwkREYDJTfZursJUf8MuvQ+0cMxu6GHYxpfhIVDuxuho9NRxdlj8MsuXox8E6ZeAYsfMuaQPG/vNI1S6gal1H6l1A9KqUke0qsrpeY7079VSsU678cqpS4opXY4jzf9taV7Qnaj2VVw3zojNLEmpPHqx8idus3hf7Yy7H//xcKOh425o+0f0i3sCeAmY0N0eGSp2jXpQ2wccEpE2iilRgD/5KIPsUMikmy2Pa2E7Eb1WnBZgtVSaGyAXz9GhSjFFmkPg/8CN/wT9i5l87yaRtqaZ2HvUh6JSIBfY6HRlWaa9utDzHn9tPN8IfCqUmWzsNTDMbtxMs2YfDyfabUkmspI9VqQPBJHYf+iaTLUieHB8MXwxjXwWhcztZjxIebKIyIO4AzQwJnWUim1XSn1lVKqp7/GdE/Ibvz8neGdr8U1RgA+jcZJ8eFZ8TSPxA2FuKF0mTSHLbecN5b3IUIptcUt13QRme52bcaHmLc8PwPNRSRTKXUV8IlSKk5EznoWUCsh+6E3sGq8UJ49ayeoA1ffYQRzHKscIuLLj4gZH2KFeTKUUhFAHeCkiAiQAyAiW5VSh4ArgC14QQ/H7EahEgrXvw8ay3D5EFNKRWL4EFtSLM8SYIzzfBiwWkREKdXQObGNUqoV0BY47Ksx/U23G66ekP7XaKzBpA+xd4APlVI/ACcxFBVAL2CyUsoB5AP3ichJX+3pb7rdcDk108MxjXWY8CGWDQz3UO5j4OPStKWVkN1IGAatekN0Q6sl0WiCglZCdqN6bePQaEIEUxPTJky4eymltimlHEqpYYEXM4Q4ugnWvVSldklrNL4IVBjoI8BdwMSKEDKkSF9rxDbv+oDVkmgqMV6dp9mQQIWBTnemFVSAjKFFvl4d05SfcvvBDiIVEQbaKzoMtAkKHIDSrjw0IYOZn1vTYaD94TQNnw4QHR1doo68vDwyMjLIzs4uS/VVgobHf6FBWDj79u4F4O1BTdjrPNdoAKKiooiJiaFataphxhGwMNCBICMjg9q1axMbG0sZN+RWfo7WgfBIOnToAEBexmk6xNS1WCiNXRARMjMzycjIoGXLquHuxcxwzIwJd0DIzs6mQYMGoauAAPo8Dg/vtFoKjU1RStGgQYMqNVrwq4Sc2/QLTbj3Ah8VmnArpQYBKKU6K6UyMCwo31JK7S6rQCGtgMBwz1lLGypqvFPV3pFAhYHejDFM05SXvUsh8wfo8WerJdFogoLeRW839q8wnJppNCGCrY1RAh2mxKOPXrtR4NA2QprQQkQsOWrWrCnF2bNnT5HrFo99WiJPeQh0faWlW7du/tMW3C3y746u+98dPeWzzujoaI/3s7KypFevXuJwOOTUqVPy2muvFUnv3bu3pKWleUzzx9133y0NGzaUuLi4IvdzcnKkZ8+ekpeXV6r6RKTUcuzcuVOaN28ur7/+epH77p+7IuQsZPny5XLFFVdI69at5bnnniuRvm/fPklKSnIdtWvXlpdeeslveX/1FlL8XSkNwHmx6L33dGglZDEFBQWSn59/8ca8VJFXr3ZdllUJvfrqqzJt2jQREUlLSyvxIhYqIU9pJ0+e9NnmV199JVu3bi1RTkTk6aefllmzZvks7wlPcvhj/fr10rVr1yL33D93RcgpIuJwOKRVq1Zy6NAhycnJkcTERNm9e7fP/I0bN5b09HSf5UtTb1VSQnpOyAOzZs3i6quvJjk5mXvvvZf8/HzS09Np374948ePJz4+nlGjRrFy5Uq6d+9O27Zt2bRpkyvPmDFjSExMZNiwYWRlZbnqrVWrFgDp6el06NCBBx54gE6dOnH06FFX2pEf09h38BBJSUnceeedrrJDhgzhqquuIi4ujunTp+OP2bNnM3jwYAAmTZrEoUOHSE5O5tFHHy2Sz1NaSkoKd9xxB6tXrzZ+qYrRq1cv6tev77HdIUOGMHv27BL3z58/z8CBA0lKSiI+Pp758+f7lcMfjRo1Yvfuogux7p+7LHKaYdOmTbRp04ZWrVoRGRnJiBEjWLx4sdf8q1atonXr1rRo0cJn+dLWW2WwSvvZtSe0Z88euemmmyQ3N1dERO6//3754IMPJC0tTcLDw2Xnzp2Sn58vnTp1krvvvlsKCgrkk08+kcGDB0taWpoAsm7dOhExhi0vvPCCq+7CXktaWpoopWTDhg1F0r7//nu5sl1bOf7zURERyczMdPWEMjMzRcQYbsTFxcmJEyeK1OlOTk6ONG7c2HVd2p6Qw+GQpUuXytChQ6V9+/byj3/8Q44dO1Ykj7eei8PhkEsvvbTE/YULF8r48eNd16dPnzZVny+GDRsmkZGRrh5G8c9dFjl79OhRZBhVeHzxxReuPAsWLJBx48a5rmfOnCkPPvigVznvvvtueeWVV/yWL029uidUhVm1ahVbt26lc+fOJCcns2rVKg4fNlzktmzZkoSEBMLCwoiLi6Nv374opUhISCA9PR2Ayy+/nO7duwOQmprKunXrPLbTokULunbtWuTe6tWrGXLrcC69zLB2cP8Vf/nll0lKSqJr164cPXqUgwcPev0MJ06coG7dsltZh4eHc9NNN7Fo0SK+/vprDh8+TPPmzdm0aZOpspGRkZw7d67I/YSEBFauXMljjz3G2rVrqVOnTpnlA1ixYoWrd1XYGyrN5/Ym59q1a9mxY0eJo1+/fq484qF36M12Jzc3lyVLljB8+EUnhN7Kl6beqoRehimGiDBmzBiee+65IvfT09OpXr266zosLMx1HRYWRuGG3OJfGm9foujoaI9tdwvbCd++BV3udd3/8ssvWblyJRs2bKBmzZr06dPHp8VsjRo1ym1Re+bMGebPn897771HtWrVeOedd0hMTDRVNicnh6ioqCL3rrjiCrZu3cqyZct4/PHH6d+/P08++aSXGnyTnZ3NX//6V5YsWcJ7773H999/z4ABA0r9uT3J2bNnzxKKCWDq1KkuRRQTE8PRoxf3dGdkZNC0aVOPbSxfvpxOnTrRuHFj1z1v5UtTb1XC1krIV5ylstbnj759+zJ48GD+/Oc/06hRI06ePOnxS+mNI0eOsGHDBrp168bcuXPp0aOH6bJ9+/blwmuTyduVR7Uu93Ly5EkgjDNnzlCvXj1q1qzJvn372Lhxo8966tWrR35+PtnZ2URFRVG7dm2vn8FTWmpqKhs2bGD48OHMnDmTtm3bmv4MmZmZNGzYsMTmyp9++on69euTmppKrVq1eP/99/3K4Y1nnnmG0aNHExsbS0JCAkuWGLuIin/ussi5du1av+137tyZgwcPkpaWRrNmzZg3bx5z5szxmHfu3LmMHDnSVPl27dqZrrdKYdU40MyckFXMmzdPkpKSJCEhQTp16iQbNmwoMbcwZswYWbBggYhcnHdIS0uTDh06yL333isJCQlyyy23yPnz511l3OeEis9TFKb9+myirLu3kSQmJsqYMWPku6OnJDs7W2644QZJSEiQYcOGSe/evWXNmjVFyhVn7NixReYxRo4cKXFxcTJx4kQRuTgn5Clt8eLFPpevR4wYIZdddplERERIs2bNZMaMGa60BQsWyCOPPFKizIoVKyQhIUGSkpIkJSVFNm/eLCIiN954o2u+qbgc7mmF7Nu3T66++mqXfPv27ZOOHS+aNLh/7rLIaZbPPvtM2rZtK61atZJnnnnGdd9d5vPnz0v9+vVLzH/5Ku/tfnGq0pyQVkIBpCyTqyV4s5fIrGGuS39L9N7Ytm2bpKamek13V0KBZOjQobJv376A12sWf5+7EKvlLC9VSQnpiWm7UZAfkHA/HTt25NprryU/Pz8AQpkjNzeXIUOG0K5du6C1WRwzn9sOcmouogzFGHyio6Pl/PnzRe7t3bvX5UcnZHmzB9RvDbd9AMDOjNMkVoA/offff58hQ4aUaxVNYx3leVeUUlkiUnJlxCJsPTEdktzneUk/0Nx1111BaUej8Ycejmk0GkuxnRKyanhoG1b8L2yfZbUUGhtT1d4RWymhqKgoMjMzq9xDLhXfLzQCIGo0HhAxfEz7s4OqTNhqTigmJoaMjAyOHz9utSiW0TY3m7Nnf+dXZ4SNX09dYO85+wau0wSfwmgbVQVbKaFq1apVmQgCZeYToX6DRtR3rnzcOOkz0qcMtFgojabiCFQs+upKqfnO9G+VUrGBFjRkKHDowIcayynPO6+Uetx5f79S6np/bflVQm6x6G8ErgRGKqWuLJZtHHBKRNoALwH/9FevxgvVaxuHRmMR5XnnnflGAHHADcDrzvq8EpBY9M7rp53nC4FXlVJKfMwwt6mdA88WiyadMBxunmacT2nhDInsxlV3wfX/AEcuPO9h2NbtIbj2ccg6CdMSSqb3fgy6/xFOH4HXu5VMv24ydB4Hv+2FGf1Kpg/8FyTdDkc3w4dDSqbfMh3aD4TDX8K8USXTR8yGVn1g32ewaELJ9Ds/gYkHSt7XaIJLmd955/15IpIDpCmlfnDWt8FbY34tppVSw4AbRGS88/pOoIuIPOSW53tnngzn9SFnnhPF6poAFL59VwFZ2IsIwOE3V3Cxo0xgT7m0TOaoCWx1u54uRoh2oHzvPIZi2igis5z33wGWi8hCb8IEKha9qXj14haLXim1RURSTLQfNLRM5rGjXFomc5iQqTzvvCld4I6ZiWkzsehdeZRSEUAd4KSJujUajf0ozztvpmwRAhWLfgkwxnk+DFjtaz5Io9HYmvK880uAEc7Vs5ZAW8Cn9a3f4ZiIOJRShbHow4F3xRmLHtgiIkuAd4APnZNQJ51C+8N/yIjgo2Uyjx3l0jKZw6dM5Xnnnfk+wpjEdgAPiohPfzKWufLQaDQasNneMY1GE3poJaTRaCylwpWQHbd8mJDpLqXUcaXUDucxPggyvauU+s1pf+EpXSmlXnbKvFMp1ckGMvVRSp1xe05li+FTOpkuV0qtUUrtVUrtVko97CFPUJ+VSZmC+qyUUlFKqU1Kqe+cMv3NQx57bLeqSAfWGJNah4BWQCTwHXBlsTwPAG86z0cA820g013Aq8F09g30AjoB33tJHwAsx7DD6Ap8awOZ+gCfBvk5NQE6Oc9rAwc8/P+C+qxMyhTUZ+X87LWc59WAb4GuxfIE9d3zdlR0T8hl/i0iuUCh+bc7g4EPnOcLgb5O828rZQo6IvI1vm2rBgMzxWAjUFcp1cRimYKOiPwsItuc5+eAvUCx/T/BfVYmZQoqzs/+u/OymvMovgoV7HfPIxWthJoBR92uMyj5z3HlEREHcAZoYLFMALc6u/ILlVKXe0gPNmblDjbdnF3+5UqpuGA27Bw+dMT4lXfHsmflQyYI8rNSSoUrpXYAvwFfiIjX5xSkd88jFa2EArblI4CYaW8pECsiicBKLv5aWEmwn5MZtgEtRCQJeAX4JFgNK6VqAR8DfxKRs8WTPRSp8GflR6agPysRyReRZAyr5auVUvHFRfZUrKLlKk5FKyE7bvnwK5OIZIqxCxjgbYzNtlZTanP4ikZEzhZ2+UVkGVBNKXVpRberlKqG8bLPFpFFHrIE/Vn5k8mqZ+Vs7zTwJYZrDXdssd2qopWQHbd8+JWp2PzBIIwxvtUsAUY7V366AmdE5GcrBVJKXVY4h6CUuhrj+5RZwW0qDGvdvSLyLy/ZgvqszMgU7GellGqolKrrPK8B9AP2Fctmj+1WQZilH4CxWnAIeMJ5bzIwyHkeBSwAfsDYY9LKBjI9B+zGWDlbA7QPgkxzgZ+BPIxfqHHAfcB9cnG14zWnzLuAFBvI9JDbc9oIXBMEmXpgDBl2AjucxwArn5VJmYL6rIBEYLtTpu+BJz18z4P+7nk69LYNjUZjKdpiWqPRWIpWQhqNxlK0EtJoNJailZBGo7EUrYQ0Go2laCWk0WgsRSshjUZjKf8fdNyzBLa0RCIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.hist(\n",
    "    x[x > cond],\n",
    "    bins,\n",
    "    density=True,\n",
    "    fill=False,\n",
    "    histtype=\"step\",\n",
    "    label=f\"empirical (t|t>{cond}) s.t. $\\lambda({cond})={lamb(cond):.3f}$\",\n",
    ")\n",
    "pl.legend(loc=\"lower left\")\n",
    "pl.xlim(xmin=0, xmax=xmax)\n",
    "pl.twinx()\n",
    "\n",
    "pl.plot(\n",
    "    bins,\n",
    "    np.where(bins < cond, 0, mixture_pdf(bins[:, None])),\n",
    "    \"--\",\n",
    "    color=\"C1\",\n",
    "    label=f\"mixture pdf\\n$\\hat\\lambda({cond})={lamb1:.3f}$\",\n",
    ")\n",
    "pl.legend(loc=\"upper right\")\n",
    "pl.ylim(ymin=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The main bias in the estimation is due to lamb'(t), independent of GMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 1) lamb(0)=1.200 true_neg_slope=1.455 gmm_neg_slope=1.366\n",
      "(1, 2) lamb(1)=0.700 true_neg_slope=1.000 gmm_neg_slope=0.908\n",
      "(2, 3) lamb(2)=0.450 true_neg_slope=0.704 gmm_neg_slope=0.774\n"
     ]
    }
   ],
   "source": [
    "for interval in [(0, 1), (1, 2), (2, 3)]:\n",
    "    true_poly = np.polyfit(\n",
    "        np.linspace(*interval),\n",
    "        [\n",
    "            np.log(lamb(b)) - sp.integrate.quad(lamb, 0, b)[0]\n",
    "            for b in np.linspace(*interval)\n",
    "        ],\n",
    "        1,\n",
    "    )\n",
    "    gmm_poly = np.polyfit(\n",
    "        np.linspace(*interval),\n",
    "        model(mx.nd.array(np.linspace(*interval))[:, None])[0].asnumpy(),\n",
    "        1,\n",
    "    )\n",
    "    print(\n",
    "        f\"{interval} \"\n",
    "        f\"lamb({interval[0]})={lamb(interval[0]):.3f} \"\n",
    "        f\"true_neg_slope={-true_poly[0]:.3f} \"\n",
    "        f\"gmm_neg_slope={-gmm_poly[0]:.3f}\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Environment (conda_mxnet_p36)",
   "language": "python",
   "name": "conda_mxnet_p36"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
